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Mass matrix 


Section area 

Correction factor for buoyancy force 
Half-beam of craft 

Crossflow drag coefficient 
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Wavelength coefficient L/A [Ca /( 2b)2 v2 


Friction drag force 


Total hydrodynamic force in x direction 
Total hydrodynamic force in z direction 
Total hydrodynamic moment about pitch axis 


Two-dimensional hydrodynamic force 
Acceleration of gravity 

Wave height, crest to trough 
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Pitch moment of inertia 
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Wave number 
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Hydrodynamic force normal to baseline 


Wave elevation r= Ty COS (kx+wt) 
Wave amplitude 


Relative fluid velocity parallel to baseline 
Relative fluid velocity normal to baseline 
Speed-to-length ratio in knots/ft!/? 
Weight of craft 


Vertical component of wave orbital velocity 
Vertical component of wave orbital acceleration 
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ABSTRACT 


A nonlinear mathematical model has been formulated of a 
craft having a constant deadrise angle, planing in regular waves, 
using a modified low-aspect-ratio or strip theory. It was 
assumed that the wavelengths would be large in comparison to 
the craft length and that the wave slopes would be small. The 
coefficients in the equations of motion were determined by a 
combination of theoretical and empirical relationships. A 
simplified version for the case of a craft or model being towed 
at constant speed was programed for computations on a digital 
computer, and the results were compared with existing experi- 
mental data. Comparison of computed pitch and heave 
motions and phase angles with corresponding experimental data 
was remarkably good. Comparison of bow and center of 
gravity vertical accelerations was fair to good. 


ADMINISTRATIVE INFORMATION 
This investigation was authorized by the Naval Sea Systems Command with initial 


funding under Task Area SR-023-0101 and completion under Task Area ZF-43-421001. 


INTRODUCTION 

Computer programs for estimating the motions of displacement ships in waves for all 
headings and speeds have been in existence for some time. Comparable computational 
schemes for planing craft do not exist except in limited and restricted cases. A program for 
planing craft would be quite useful to the small craft designer, providing a means for 
systematically exploring the effects of numerous design variations on performance of the 
craft in waves. With minor modification, the program could also be used to examine the 
merits of a hybrid craft design, e.g., a combination of planing craft and hydrofoil. 

Predicting the motions of a planing craft in wave’s is by no means a simple problem. 
The analytical description of a high-speed craft, planing in waves, involves several different 
types of flow phenomena, including planing; hydrodynamic impact, and, to a lesser extent, 
surface wave generation and hydrostatics. Also, the mathematics tend to become nonlinear 
rapidly as the motion increases or, like the real craft, can in some instances exhibit large 
instabilities such as porpoising. 

Development of a computer program that would take into account all of the previously 
described factors and would be applicable for a wide range of speed and wave conditions 
requires a careful and systematic study in several stages with appropriate verification at each 


stage. To lay the foundation for such a general program, a simpler problem has been 


formulated in this report with potential for expansion and generalization to the more 
complicated case. The simpler problem is that of a V-shaped prismatic body with hard chines 
and constant deadrise planing at high speed in regular head waves. 

The mathematical formulation is analogous to low-aspect-ratio wing theory with 
provisions for including hydrodynamic impact loads, essentially a strip theory. Surface wave 
generation and forces associated with unsteady circulatory flow are neglected, and the flow is 
treated as quasi-steady. The mathematical formulation is an empirical synthesis of several 
theoretically derived flows describing the overall craft hydrodynamics. Wave input is restricted 


to monochromatic linear deepwater waves with moderate wavelengths and low wave slopes. 


MATHEMATICAL FORMULATION 
GENERAL 
Consider a fixed coordinate system (x,z) (Figure 1) with x axis in the undisturbed free 
surface, pointing in the direction of craft travel, and the z axis, pointing downward. If the 
motions of the craft are restricted to pitch 0, heave Zcq> and surge Xcq, the equation of 


motions can be written as 


MXcg = T,,-N sin @.—Djcos 0 
MZcc = T,-Ncos#+Dsin@ +W 


10 = Nx, + Dx Tx, (1) 


where M is mass of craft 
I is pitch moment of inertia of craft 
N_ is hydrodynamic normal force 
D_ is friction drag 
W is weight of craft 
is thrust component in x direction 
is thrust component in z direction 
x, is distance from center of gravity (CG) to center of pressure for normal force 
Xq is distance from CG to center of action for friction drag force 
x. is moment arm of thrust about CG. 


Equation (1) is exact; however, defining the hydrodynamic forces and moments in waves 


can be extremely difficult. 


A high-speed craft moving in waves may transit through several regimes that have 
different hydrodynamic flow characteristics. For example, as the craft moves away from the 
crest of wave, the flow may be characterized by unsteady-state planing until the craft collides 
with the oncoming wave crest and enters another regime in which impact forces are important. 
After the impact, the craft may enter still another regime in which it is planing but in which 
buoyancy forces are rather significant. 

The most promising approach to a method that would incorporate all three types of flow 
conditions into a general formulation would seem to be a modified strip theory. The 
mathematical justification for this approach is not rigorous; however, there is sufficient 
precedent to expect promising results. For example, impact loads on landing seaplanes can 
be estimated reasonably well using a strip theory incorporating the Wagner two-dimensional 
(2-D), expanding-wedge theory,! and Chuang? has provided a strip method for determining loads 
on an impacting prismatic form that agrees extremely well with experimental results. 

More recently, Martin? has developed a linear strip theory for estimating motions of a 
planing craft at high speed, which shows good agreement with experimental results. A 
nonlinear model of the equations of motion would be expected to provide, in addition to the 
motions, reasonable estimates of the vertical accelerations which are an important consideration 


in designing a planing craft. 


TWO-DIMENSIONAL HYDRODYNAMIC FORCE 

Implicit with any strip method is the need to define the 2-D hydrodynamic force acting 
on an arbitrary cross section of the body. The 2-D flow problem is not simple; however, it 
lends itself to an empirical approach, using a combination of techniques used in hydrodynamic 
impact and low-aspect-ratio theories. 

The typical cross section of a hard-chine, V-shaped prismatic body such as that being 
considered here is shown in Figure 2. Figure 2 actually illustrates two different idealized- 
flow conditions, assumed to represent the crossflow during unsteady planing, depending upon 
whether the flow separates from the chine (Figure 2a) or not (Figure 2b). Nonwetted-chine 
flow conditions are typical of the sections near the leading edge of the wetted length of the 
craft. Wetted-chine flow conditions are more typical of sections near the stern, except 
possibly in the most extreme motion and wave conditions. Some sections between leading 
edge and stern may alternate between flow conditions as the wetted length changes with the 


motions. 


*A complete listing of references is given on page 33. 


3 


The normal hydrodynamic force per unit length f, acting at a section, is treated as 
quasi-steady and is assumed to contain components proportional to the rate of change of 


momentum and the velocity squared (drag term), i.e. 
Pes all Gate: 2 
=~ |p (ma) + Ope bVv (2) 


where V is the velocity in plane of the cross section normal to the baseline 


m, is the added mass associated with the section form 


C,, . is the crossflow drag coefficient 


D,c 
p is the density of the fluid 
b is the half beam. 
For sections near the leading edge of the wetted length with nonwetted chine, the 
added mass is assumed to be defined in the same manner as during an impact which for a 


V-shaped wedge is given by 
m, =k, 7/2 pb? (3) 


where k, is an added-mass coefficient that may also include a correction for water pileup- 
k, is assumed to be 1.0 without pileup correction. 


The rate of change of momentum of the fluid at a section is given by 


io) 
dg 


dé 


D alee 
ay DS Pid eae (Ca De (4) 


where & is the body coordinate parallel to the baseline; see Figure 1. The last term on the 
right-hand side of Equation (4) takes into account the variation of the section added mass 
along the hull. This contribution can be visualized by considering the 2-D flow plane as a 
substantive surface moving past the body with velocity U = -dé/dt tangent to the baseline. 
As the surface moves past the body, the section geometry in the moving surface may change 
with a resultant change in added mass. This term exists even in steady-state conditions and 
is the lift-producing factor in low-aspect-ratio theory. 

The added mass of a section with fully wetted chines has not been developed to the 


same extent as the V wedge. In steady-state planing problems such as those of Shuford,* 


the crossflow is treated as a Helmholtz-type flow in which the Bobyleff results are used for 
estimating drag coefficients. Helmholtz flows are applicable only to steady-state conditions; 
so, it is assumed that the added mass for the fully wetted chine flow can be determined from 
Equation (3) using the value of the half-beam at the chine. In using the Shuford approach, 
it is assumed that the crossflow drag coefficient for a V-section is equal to the drag of a flat 


plate (C, , = 1.0) corrected by the Bobyleff flow coefficient approximated by cos 8, i.e. 


Coc = 1.0 cos 8 (5) 


The Bobyleff flow coefficient is the theoretical ratio of the pressure on a V-section to that 
experienced by a flat plate for a Helmholtz-type flow. 

The same approximation is used for estimating the drag coefficient for nonwetted chine 
sections, using the instantaneous value of the half-beam at the free surface. 

An additional force acting on the body is the buoyancy force f,,- This force is assumed 
herein to act in the vertical direction and to be equal to the equivalent static buoyancy force 


multiplied by a correction factor, i.e. 
fi =tee(A) (6) 


where A is the cross-sectional area of the section, and a is a correction factor. 

The full amount of the static buoyancy is not realized because at planing speeds the water 
separates from the transom and chines, reducing the pressure at these locations to atmospheric 
or less than the equivalent hydrostatic pressure. A greater reduction is realized in the 
buoyancy moment because of the corresponding shift in the center of pressure. Shuford4 

in his work on steady-state planing recommended a factor of one-half to obtain the correct 
buoyancy force. In the following computations, the buoyancy force was corrected by a 
factor of one-half, ie., a= 1/2. The buoyancy moment, computed as the static buoyancy 
force multiplied by its corresponding moment arm, was corrected by an additional factor of 
one-half to obtain the proper mean-trim angles. 

Equation (2) is a synthesis of several idealized flow conditions combined in an empirical 
manner. In all of these flows, it is assumed that the net relative movement of the fluid past 
the body is in an upward direction. This condition may not always be met in the case of 
unsteady planing in waves. Closer scrutiny will be required to determine what limitations 
will be imposed upon the problem as formulated and/or what modifications will be required 


to improve the formulation. 


TOTAL HYDRODYNAMIC FORCE AND MOMENT 

The total normal hydrodynamic force acting on the body is obtained by integrating the 
stripwise, 2-D, hydrodynamic force given by Equations (2) and (6) over the wetted length 2 
of the body. A body coordinate system (£,¢) with its origin at CG and the & axis pointing 
forward parallel to the baseline of the body is defined in Figure | to facilitate this integration. 
The hydrodynamic force acting in the vertical or z direction of the fixed integral coordinate 


system is given by 


-Ncos@) = FU(t) = / f cos @dé +f igat 
Q Q 


= | ftraeoveo + m, (,t) V(E, t) 
Q 


a0} 
dg 


+Cy Deb EN VE,D} cos dé 


+apsade] (7) 


where the integration is taken over the instantaneous wetted length. Similarly the force 1 


acting in the horizontal or x direction is given by 


F, = [tin oat 


= — [ {mge.0 ven + m,(é,t) V(t) 


- UE) = [m,(E.0 VED) 


+ Cp .é.teb,0V7(E,b} sin 0 dé (8) 


Wave forces are obtained by neglecting diffraction and assuming that the wave excitation 
is caused both by the geometrical properties of the wave, altering the wetted length and 
draft of the craft, and by the vertical component of the wave orbital velocity at the surface 


w,, altering the normal velocity V. The horizontal component of orbital velocity is neglected, 


since it is assumed small in comparison with the forward speed arae The velocities U and V 


may then be written as 


U= Xog Cos 6 - (Z¢g-Wz) sin 0 


V = Xcq sin @ - bg + (Z6g-Wz) cos 9 


The depth of submergence h of the body at any point P(£,¢) may be determined by 


h = Zceg - & sin @ + § cos @ -r 


where r is the instantaneous value of the wave elevation directly above the point. 


For regular head waves the wave elevation for a linear deepwater wave is 


TS atcOs k(x+ct) 


where r, is the wave amplitude 


k is the wave number 
c is the wave celerity. 
At point P(é,¢) 


X = Xo, t Ecos 8 +§ sin 8 


(9) 


(10) 


(11) 


(12) 


The hydrodynamic moment Fg about CG is obtained in a similar manner by integrating 


over the wetted length the product of the normal force per unit length and the corresponding 


moment arm. 


Fees - [re nece =i cos 6 &dé 
Q Q 


=f {me001E0 + ieOVED 


— ULE 3 (mE OVE) + Gy EDP bE DVED 


+apgA cos 6} Edé (13) 


EQUATIONS OF MOTION, GENERAL 
Integrating the first term in Equations (7), (8), and (13) provides hydrodynamic forces 
and moments proportional to acceleration of the motion. These can be combined with the 


inertial terms of the rigid body to give the following equation of motion 


(M+M, sin? 9) Xa +(M, sin @ cos Nae =(QFysin 0)0 
= 1, +F, -Dcosd (14) 


(M, sin 6 cos 9) Xa +(M+M, cos” Ware -(Q, cos 0)0 


=17,+F,+Dsin@+W 


-(Q, sin 0)X og -(Q, C08 M%og + (1+1,)4 


= Fg -Dxg + Tx, 


where M,(t) = [ me.nay 
Q 


Q,(t) = of m,(£,t)Edé 
Q 


L(t) = We m, (£,t) £7 dé 
Q 


F, = F, -{-(M, sin? )%¢q -(M, sin @ cos 0)%og + (Q, sin 0)6} 
Ee alse {appropriate acceleration terms} 
Fg = Fo - {appropriate acceleration terms}. 


A detailed evaluation of the integral expressions for the hydrodynamic forces and 


moments is provided in Appendix A. 


The solution to Equation (14) is cumbersome; however, it can be accomplished using 


standard numerical techniques. Introducing the state vector [x;, X4,X3,X4. Xe, X¢] 


where X) = Voc 
Xy = 20g 
or 
ACG 
CELA 
X¢ = 


Equation (14) can be rewritten, using matrix algebra, as 


> > 

Ax = g (15) 
so that 

Te eA eal ee 

x=A’g (16) 


where A“! is inverse of the inertial matrix A. Equation (16) is now in a form that lends 
itself to integration by using a numerical method such as the Runge-Kutta-Merson integration 


routine. 


EQUATIONS OF MOTION, SIMPLIFIED FOR CONSTANT SPEED 
Assuming that the perturbation velocities in the forward direction are small in comparison 
to the speed of the craft, the equations of motion may be further simplified by neglecting 


the perturbations and setting the forward velocity equal to a constant, i.e. 


og = CONSTANT 


If it is also assumed that the thrust and drag forces are small in comparison to the hydrody- 
namic forces and that they are acting through the center of gravity, the equations of motion 


may be written as 


Kats iO 
(M+M, cos? 0)%., -(Q, cos0)0 = F, +W 
-(Q, C08 O)Zcg + (I+1,)6 = Fo 


These equations also represent the case of the craft (model) being towed through CG at 
CONSTANT speed. Based upon the previously described equations of motion, a computer 
program has been written in FORTRAN language to compute the motions of a prismatic body, 
planing in regular head waves at high speed. A listing of the program along with the 
appropriate flow chart is presented in Appendix B. The listing contains reference to thrust 
and drag terms; however, they have no significance, except to provide a starting point for 


possible updating of the program to include these terms in the future. 


COMPARISON OF COMPUTED RESULTS WITH EXPERIMENTS 

Computations of pitch and heave motions and heave and bow accelerations were made, 
using the computer program for comparison with the experimental results of Fridsma.> 
Fridsma tested a series of constant-deadrise models of various lengths in regular waves to 
define the effects of deadrise, trim, loading, speed, length-to-beam ratio and wave proportions 
on the added resistance, heave and pitch motions, and impact accelerations at the bow and 
center of gravity. Figure 3 shows the lines of the prismatic models. The models were towed 
at CG with a system that permitted freedom in surge. The computer program simulates the 
model being towed at constant speed with CG at the baseline. 

Table 1 presents some characteristics of the model and experimental conditions for 
which comparisons were made. Most of the comparisons have been made at a speed-to-length 
ratio Wile of 6.0 where the mathematical model is expected to be most representative. A 
limited comparison has also been made at Wifes = 4.0; however, no comparison has been 
made at Wilx/ ale = 2.0. At this speed, the model (or craft) operates in the displacement mode 
for which the mathematical formulation is not valid. 

The average computer run corresponded to 10-second, real-time, model scale; however, 
only the last 2 seconds were considered free of transient effects. An example of the compu- 
ter time histories of pitch and heave motions is shown in Figure 4. Although the motions 
are periodic, they are not perfectly sinusoidal; consequently, in determining phase relationship, 
the peak, positive-pitch value (bow up) and the peak, negative-heave value (maximum upward 
position of CG) were used as reference points. There was a difference when the opposite 


peaks were used. 


TABLE 1 — MODEL CHARACTERISTICS AND WAVE 
CONDITIONS FOR COMPUTATIONS 


(Model Length = 114.3 cm (3.75 ft); L/b = 5; Ca = 0.608) 


CONFIGURATIONS 


B LCG Radius of 
SYMBOL Gyration V/V L 
40 


20 
20 
10 
30 


WAVE CONDITIONS FOR CONFIGURATION —-— 


ee | eee 


Corresponding time histories of bow and CG accelerations are shown in Figure 5. The 
bow acceleration was computed at Station 0. As can be seen in these plots, the impact 
accelerations ranged in magnitude from cycle to cycle. The maximum impact (or negative 
value) acceleration computed during the final 2 seconds of run was used in the comparisons 
with experimental values. In some instances, particularly near resonance, the maximum 
impact acceleration was more than twice the average impact value. 

Figure 6 shows a comparison of variation of computed and experimental pitch and heave 
motion with wave height for the 20-degree deadrise model in a 15-foot wavelength and for a 
speed-to-length ratio of 6.0. Figure 7 shows the corresponding impact acceleration at the 
bow and CG. The computed results closely follow the experimental data, except for CG 
acceleration at the extreme wave height condition, where the computed value is apparently 
much lower. Experimental data show that the model was leaving the water at this wave- 


height condition. The computer model did not leave the water but came very close; 
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see Figure 8. Figure 8 is a trajectory of the computer model relative to the wave for a 
selected cycle of motion. The computer model behaves very much as expected. On the left- 
hand side of the figure, the craft is planing down the crest of the wave and, as it approaches 
the wave trough, comes very close to leaving the water before slamming and submerging 
itself deeply into the front of the oncoming wave crest. 

Figures 9 through 14 show comparisons of the computed and experimental pitch and 
heave motions at VIS L = 6.0 through a range of wavelengths and at a constant wave height 
of 2.54 centimeters (1 inch) for deadrise models with 10, 20, and 30 degrees. The data have 
been plotted with respect to the coefficient C), defined by Fridsma as L/A [Ca/(L/2b)7] M3) 
Note that in our notation, b is the half-beam. 

Comparisons of heave and pitch for the 10-degree deadrise model shown in Figures 9 
and 10, respectively, show excellent results. The computer model accurately predicts the 
secondary peaks in the pitch and heave responses at C, = 0.19. At this condition, the physical 
experimental model rebounds so as to fly over alternate waves. The computer model oscillates 
at half the wave-encounter frequency and comes close to leaving the water at alternate 
encounters with the wave. It does not quite leave the water to fly over alternate wave crests; 
nonetheless, it is a good representation of the actual motion. 

The heave and pitch comparison for the 20-degree deadrise model at WH TE = 6.0 is also 
excellent as can be seen in Figures 11 and 12, respectively. No experimental phase data for 
the condition were reported for C) greater than 0.072; however, extrapolated results (not 
shown) are in line with the computed results. The pitch and heave results shown in Figures 
13 and 14 for the 30-degree deadrise model are good; however, responses at C) = 0.048 and 
C) = 0.072 are higher than the experimental results. 

For practical considerations a computational scheme for planing boat motions should be 
valid for a range from approximately Wale = 4.0 to Vis = 6.0. Computations of the 
motions were made for W/E = 4.0 for the 20-degree deadrise model; see Figures 15 and 16. 
Again the comparison of the computed heave and pitch response with experimental results is 
excellent. 

Comparisons of the computed and experimental impact accelerations (or largest negative 
values) are presented in Figures 17 through 20. Figures 17 and 18 show bow and CG 
accelerations for the 10-degree deadrise model; Figure 19 shows similar results for the 20- 
degree deadrise model; Figure 20 shows the results for the 30-degree deadrise model. In all 
cases, the comparison appears to be fair to good. In the shorter wavelengths, A/L = 1.0 and 
A/L = 1.5, the computed accelerations are higher than the corresponding experimental values. 
This is most pronounced for the 10-degree deadrise angle model. 
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CONCLUSIONS AND RECOMMENDATIONS 

A mathematical model of a craft having a constant deadrise angle, planing in regular 
waves, has been formulated using a modified low-aspect-ratio or strip theory. It was assumed 
that the wavelengths were long in comparison to the craft length and that the wave slopes 
were small. The coefficients in the equations of motion were determined by a combination 
of theoretical and empirical relationships. 

A simplified version for the case of a craft or model being towed at constant speed was 
programed for computations on a digital computer, and the results were compared with 
existing experimental data. 

The comparison of the computed pitch and heave motions and phase angles with the 
corresponding experimental data gave remarkably satisfying results. Comparison of the bow 
and CG accelerations was fair to good. 

In summary, the previously described mathematical model appears to be a valid represen- 
tation of a planing craft in waves for the specific craft geometry and wave conditions 
considered. 

To make the computer program more valuable to the designer the following additional 
work is recommended: 

1. Improve estimates of hydrodynamic coefficients to obtain better acceleration data 
and to include more complicated ship geometry. 

2. Determine added resistance in waves. 

3. Include freedom to surge and to add components of propulsion. 


4. Extend to the case of irregular waves. 
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Figure 1 — Coordinate System 


Figure 2a — Flow Separation from Chine 


Figure 2b — Nonwetted Chine 


Figure 2 — Types of Two-Dimensional Flow 
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Figure 4 — Sample Time Histories of Computed Pitch and Heave Motions 
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Figure 5 — Sample Time Histories of Computed Accelerations 
of Bow and Center of Gravity 


17 


DOUBLE AMPLITUDE OF HEAVE (cm) 


DOUBLE AMPLITUDE 
PITCH MOTION (deg) 


(inches) 


LEAVING WATER 


A EXPERIMENTAL (REFERENCE 5) 
O COMPUTED 


0.1 0.2 0.3 
WAVE HEIGHT/BEAM 


Figure 6 — Variation of Pitch and Heave with Wave Height 
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Figure 7 — Variation of Acceleration of Bow 
and Center of Gravity with Wave Height 
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Figure 9 — Heave Response for 10-Degree Deadrise Model at V/,/ L= 6.0 
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Figure 10 — Pitch Response for 10-Degree Deadrise Model at V/,/ L = 6.0 
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Figure 11 — Heave Response for 20-Degree Deadrise Model at V//L= 6.0 
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Figure 12 — Pitch Response for 20-Degree Deadrise Model at V/,/ L = 6.0 
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Figure 13 — Heave Response for 30-Degree Deadrise Model at V/,/ L = 6.0 
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Figure 14 — Pitch Response for 30-Degree Deadrise Model at V/,/ L= 6.0 
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Figure 15 — Heave Response for 20-Degree Deadrise Model at V/,/ L = 4.0 


PITCH (0_/27H/)) 


p 


A EXPERIMENTAL (REFERENCE 5) 
O COMPUTED 


0 0.05 0.10 0.15 0.20 0.25 0.30 
Cy 


Figure 16 — Pitch Response for 20-Degree Deadrise Model at V/,/ L= 4.0 
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Figure 17 — Bow Acceleration for 10-Degree Deadrise Model at V/,/ L = 6.0 
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Figure 18 — Center of Gravity Acceleration for 10-Degree 
Deadrise Model at V/,/ L = 6.0 
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Figure 19 — Bow and Center of Gravity Accelerations for 20-Degree 
Deadrise Model at V/,/ L = 4.0 and V/,/ L= 6.0 
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Figure 20 — Bow and Center of Gravity Accelerations 
for 30-Degree Deadrise Model at V/,/ L= 6.0 
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APPENDIX A 
EVALUATION OF HYDRODYNAMIC FORCE AND MOMENT INTEGRALS 


The hydrodynamic force the craft experiences in the vertical direction as derived in the 


text is: 


‘ dm,V ; ; 
F,=-f {mV-u 9E MV iste IDV }cosoat+ J apeadt 


where U = XGG cos @ - (z-w,) sin 0 
and 
V = Xcg sin ®@ + (z-w,) cos 8 -6§ 


Another force acting in the vertical direction is the weight of the craft. 
The first two terms of the integral are evaluated by making the substitutions 


V =Xeq sin? —-@€ + Z—, cos'8-w, cos @ 


cE) ars cos 6 - ZeG sin 0) +w, 6 sin 0 


. Ow 

a = -8 - = cos 8 
A ONE) ; 
= aa 
o£ 0& 
dw ow 

aah z 
Taha he ouRE 


and noting that 


dUV 
= oe 
on Jom cine 


Using the previously described substitutions, the force becomes 


am, 
fov BE d§=-UVm, 
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cos 0) 


aoe {-M, cos 8 Z., - M, sin 6 nee + Q,6 + M, 0 (Zeg sin 6 - ets 


dw, A 

+ fm, a cos 0 ag- ['mw,0 sin 6 dé 

Q Q 

dw, vi dw, 
= V — si + —_ 
fm aE sin 6 dé A m,U DE cos 6 dé 
-UV - | Vm,dé - 2 
mM, fight i m, dé pf ca,cb¥ dé} cos 6 

+ [ apenct 

Q 


where M, = [mae 
Q 
and 
Q, = [meat 
Q 


This is essentially the form in which the integrals have been computed in the program. 
The rate of change of the sectional added mass in the third term of the integral 
expression is derived by relating it to the rate of change of depth of fluid penetration of the 


section. The added mass of a section is assumed to be equal to 
ma aka me pb? 

for which the time derivative is 
m, = k, tpbb 


where b is the instantaneous half-beam of the section, and k, is an added-mass coefficient, 
assumed to be constant. A value of k, = 1.0 was used in the computations contained in this 
report. For sections with constant deadrise, which is an imposed limitation of this work, 


the half-beam is related to the depth of penetration by 


b=dcotB 
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where d is depth of penetration, and @ is deadrise angle. 
Taking into account the effect of water pileup, the effective depth of penetration d, is, 


according to Wagner 


d, =7/2d 


e 


and 


b=d, cot 6 = m/2 d cot B 


where 77/2 is the factor by which the wedge immersion is increased by the pileup. Using this 


expression for the half-beam, the rate of change of sectional added mass becomes 


m, = kampb(n/2 cot B)d 


This expression is valid for penetration of the section up to the chine. When the immersion 


exceeds the chine, the sectional added mass is assumed to be constant, i.e., 


2 


may = koma plOme 


rh, = 0 


where Dee is the half-beam at chine. 


The submergence of a section in terms of the motions is given by 
h=z-r 


where z = Zo, - & sin 0+¢€ cos 0 


ia COS (kGee +£ cos @ +f sin 9) + wt} 


For wavelengths which are long in comparison to the draft and for small wave slopes, the 


immersion of a section measured perpendicular to the baseline is approximately 


Dial 


~ 


cos @ -ysin#@ 


where v = wave slope 
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The rate change of submergence d is given by 


z-Tf + (z -1r) d(cos 8 - usin 0) 


eo" ee 


cos 6 - usin 0 (cos @ - v sin 6)2 ot 


Since immersion (z-r) is always small in the valid range of the previously described 


expression, the relationship can be further simplified to 


WG gee 


~ cos 6 -vsin 8 
and 
(z -1) 


m, ~ k,tpb(n/2 cot 8) ~—— 


The expansion of the integral expression for the hydrodynamic moment in pitch 


follows the procedure used for the vertical force. The results are summarized as follows 


Fo = -1,0 + Q, cos O Zog - Q,4Z¢g sin 8 - kg cos 8) 


dw, z 
~ fm, cose =? kat + [m6 sin 0 w, dé 
Q ! Q 
+ [ vingeat + [ocybveede 
Q Q 
+m,UVé| + [m,vuee 
stern Q 
uf dw, 
st; m,V —— sin @éd 
us dw, 
- / m,U = cos0@éd 
J m,U Gp cos 0 Edt 


+ [aoe cos 6 €dé 
Q 


The only additional moments are the buoyancy moments. All other moments are considered 


to be zero for the specific problem considered in this report. 
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APPENDIX B 
COMPUTER PROGRAM DESCRIPTIONS 

OVERVIEW 

The equations of motions developed in the previous sections of this report have been 
solved by means of digital computer programs. Two major programs have been developed: 
the first (MAIN) solves the equations of motion using the Runge-Kutta-Merson integration 
algorithm and generates time histories that are stored on the system disk. The second 
(PLTHSP) generates California Computer Products Company (CALCOMP) pen plots from the 
disk files. All programs were designed to operate on the Control Data Corporation computer 
system, located at the David W. Taylor Naval Ship Research and Development Center in 
Carderock, Md. 

Descriptions of input data required to execute the programs, job control cards, and 
programs follow. Sufficient detail is presented for this appendix to serve as a manual for 


use and maintenance. 


JOB CONTROL CARDS FOR PROGRAM MAIN 
Job control cards for program MAIN which computes time histories of the motion 
variables, are described as follows. If CALCOMP plots are not desired, TAPES need not be 
cataloged. 
Job Control Language Card: Comment 
Job Card 


Charge Card 


Standard facility card 
Standard facility card 


REQUEST,TAPE9,*PF. 
REQUEST,TAPE2,*PF. 
REQUEST,TAPE4,*PF. 


ATTACH,BINAR,SEFZARNICKNEWB, 
ID=XXXX. 


ATTACH,NSRDC. 
LDSET(LIB=NSRDC). 
BINAR. 


REWIND,TAPE2. 
REWIND,TAPE4. 


COPY(TAPE2,OUTPUT) 
COPY(TAPE4,OUTPUT) 


Bg, 


Reserves space for CALCOMP plot data 
Print output file 1 request 
Print output file 2 request 


Attaches binary run file 


Attaches library routines 
Loads library routines 


Loads and executes run file 
Rewinds time-history files for printing 


Prints time-history file 


Prints time-history file 


Job Control Language Card: Comment 


CATALOG,TAPE9, SEFZARNICKDATA.., Catalogues file for plot. 
ID=XXXX. (SEFZARNICKDATA CAN BE ANY NAME) 


7/8/9 END OF RECORD 
DATA CARDS (1-5) 
6/7/8/9 END OF FILE 


INPUT DATA CARDS FOR PROGRAM MAIN 

Input data used by program MAIN are read from data cards in NAMELIST and in 
standard format. A description of the FORTRAN symbols appearing in NAMELIST follows. 
For simplicity in the text that follows, it is assumed that NAMELIST input occupies only 


one card. More cards can be used if necessary. 


Card 1(NAMELIST FORMAT, / i) 
A The absolute error for KUTMER (six values) 
NPRINT If=1, print normal output 


If=2, matrix, inverse matrix, F-column matrix, and KUTMER results 
If=3, integral results 


If=4, calculated values constant for given input values 


NPLOT If=0, no plot 
If=1, printer plot of results 
END Number of runs to be made 
W Weight of craft in pounds 
BL Boat length in feet 
TZ Thrust component in z direction 
TX Thrust component in x direction 
XECG Distance from center of gravity to center of pressure for drag force in feet 
XP Moment arm of propeller thrust 
XD Distance from center of gravity to center 
DRAG Friction for drag force 
RO Wave height 
LAMBDA Wavelength 
RG Radius of gyration in feet 
a Propeller thrust in pounds 
GAMMA Propeller thrust angle in degrees 
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Card 1 (continued) 


ECG Longitudinal center of gravity 
NCG Vertical center of gravity, nondimensionalized by ship length 
KAR Added-mass coefficient 
BETA() Dead-rise angle in degrees 
EST() Station position in feet 

NUM Number of stations 

XA Initial time 

XE Stop time 

HMIN Minimum step size 

HMAX Maximum step size 

EPS Error criterion 


Card 2 (Format 8F10.0) 
(X(1),I=1,6) Initial conditions 


X(1) Velocity 
X(2) Z 

X(3) i) 

X(4) xX 

X(5) Z 

X(6) 0 deyrees 


Card 3 (8F10.0) 
START Time to turn on (RMP) function (see page 48) 
RISE Duration of RMP 


Card 4 (8F10.0) 


TME Time at which integration interval is to be changed* 
HMX New maximum interval size after TME 
HMN New minimum interval size for KUTMER to subdivide 


*If this option is not used set TME to stop time on run. 


4] 


Card 5 (8F10.0) 


PERCNT Percentage of boat length subtracted from longitudinal center of gravity to 


obtain X - point where acceleration computations are made 


JOB CONTROL CARDS FOR PROGRAM PLTHSP 
Job control cards for program PLTHSP which generates CALCOMP plots of time 


histories computed by program MAIN are described in this section. 


Job Control Language Card: 
Job Card 
Charge Card 
REQUEST,TAPE7,HI. 


Comment 
Standard facility card 
Standard facility card 
Tape for CALCOMP plot data 


VSN(TAPE7=CK0323). Volume serial number of tape for 
CALCOMP plot 

ATTACH,CALC936. Attaches CALCOMP library routine 

ATTACH, BINAR,SEFZARNICKPLOTB, Attaches plot program run file 

ID=XXXX. 

LDSET(LIB=CALC936) Loads CALCOMP library routines 

BINAR. Runs plot program 

7/8/9 END OF RECORD 

DATA CARDS 


6/7/8/9 END OF FILE 
INPUT DATA CARDS FOR PROGRAM PLTHSP 
Two or three data cards are made ready by PLTHSP, depending on the options selected. 


Standard input format is employed. A description of the necessary data cards follows. 


Card 1 (8F10.0 Format) 


XAXIS Length of x axis in inches 

YAXISP Height of pitch component axis in inches 
YAXISH Height of heave component axis in inches 
latit Height of lettering in inches 


Card 2 (110 Format) 
IA If=0, no plots for bow acceleration and center of gravity acceleration 


If=1, plots previously mentioned information 
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Card 3 (8F10.0 Format) - Only Necessary If |A=1. 
YAXISB Height of bow acceleration axis in inches 
YAXISC Height of CG acceleration axis in inches 


PROGRAM MAIN 

Program MAIN reads all necessary input data from cards, sets up initial values, computes 
constants, calls KUTMER to determine the state variables at TIME for the period from XA 
to XE in increments of HMAX. A table state variables is created for every PTIME-th value. 
The values for \/H and 0/27 H/d are calculated and printed. If the plot option is on, a 


printer plot will be produced. 


Subroutine COMPUT(X) 
This routine computes pitch moment NL and lift force FL, excluding added mass terms, 
using values of integrals computed in subroutine FUNCT. The argument X contains the state 


vector. 


Subroutine DAUX 
This subroutine is called from KUTMER or EULER. It determines the values of m,. b, 


and bl*, based on the following equations 


hed) = Zeg = E(I) sin 8 + (1) cos 6 - r(1) 


where r(I) = r, cos k ixee + €(1) cos 0 + €(I) sin @ + ct] 


Then for 


hy (1) 0) 


hy, 


OD) = ae at) Ga 


where V(I) = -rk sin 0 ixee + &(1) cos 6 + (1) sin 6 + ct] 


If 
d(1) > b,,() tan (BC) 2/n) 
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set 


m, (1) * Mamax 1) 
b(1) = b,, WD 
b1(1) =0 


Mamax(l) = k(1) (p/2)mb? (1) 


If 

d(1) <b, (D) tan (8(1)) (2/7) 
set 

b(D = d(D) cot (B(1)) (7/2) 

b1d) = bd) 

m, (1) = k,(1) (p/2)rb?(1) 
for 


hy, () <0; 
ma) ES OF b() =0, bl(l)=0 


This subroutine then calls FUNCT which in turn calls COMPUT to determine the values 
of N, and F the lift force and moment. The values of Ny and F are used to compute 


the following 


F, = T, +F, sin@ -D cos 0 
Elis Beet Ecost0kn Disinig +W 
F, 


=N, -Dyg + Tx, 


*b1 array is set up for integrations for portion of hull for which chine is not immersed. 
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The mass inertia matrix is 


A, =M+M, sin? 0 


Ap = M, sin 6 cos 8 


A,, = -Q, sin 6 
Ay, = Ay 
Ago =M +M, cos” 6 


A3, = Ay 
Aj. = Ag; 
A,, ~I1tI, 


The matrix is inverted by the system routine MATINS. The inverted matrix is then 


used to solve the following equations which determine the state vectors. 


om ae ee i “il 
Xog =A Fy t Ay” F, + Ay3” F; 


ed eM emer el =i Z 
Zog = Agy F, t+Agy F, +A; F; 


gerbils OE “1 “i 
@ =A,, F, +A,. F,+A;, F, 


Subroutine FUNCT (X) 
This routine evaluates various integrals appearing in the force and moment mathematical 
models. The integrals are evaluated, using a trapezoidal integration algorithm. The argument 


x contains the state vector. A list of integrals that are evaluated is presented. 
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Subroutine INPUT 

This routine reads in NAMELIST/HSP/ which contains the initial data concerning the 
craft and sea conditions pertinent to all the runs to be made. It is set up so that most of 
the data are given default values by means of data statements in subroutine INPUT. These 
data statements can be overridden during execution by reading values in on cards. For further 
explanation of the specific variables see section on the input data cards. 

This routine also “‘initializes’’ constant such as 7, p, and g. It uses the input values to 
calculate the keel profile and planform arrays, NO and BM, wave constants, system mass and 


inertia, and maximum mass and depth of chine at each station. 
Subroutine KUTMER (NEOS, TIME, HMAX, X, EPSE, A, HMIN, FIRST) 


This is a Runge-Kutta-Merson integration routine that is capable of changing the size of 


the interval over which it integrates to meet specified error criteria. It is therefore an 


46 


accurate method for a system that may oscillate more rapidly than the initial integration 
interval. A minimum step size prevents the routine from subdividing the interval indefinitely. 


The input arguments are: 


NEQS Number of dependent variables in the x array 

TIME Actual time (independent variable) 

HMAX Increment for which the solution is to be returned 

X Vector of dependent variables 

EPSE Relative error criteria specified for each component of x and used for the 
components of x less than the absolute value of A 

A Absolute error criteria 

HMIN Minimum step size allowed 

FIRST Set to zero on first call; a value of | is assigned by KUTMER on subsequent 
calls for which the error criteria are satisfied, otherwise a value of 2 is 
assigned 


Subroutine PLOT2 (F, FMIN, FMAX, NVAR, NFUN, N1, N, XO, DELX) 
Data stored in the two-dimensional array F are plotted, using the printer by subroutine 
PLOT2. As many as 26 different functions, having evenly spaced abscissa values, can be 


plotted. The output is written on Unit 6. A description of variables follows. 


F Array containing data to be plotted; the Jth point of the Ith function is 
stored in F(I,J) 

FMIN An array of minimum functional values; the minimum of the Ith function 
is stored in FMIN(1) 

FMAX Same as FMIN only for maximum values 

NVAR An array of titles for each function to be plotted 

NFUN Number of functions to be plotted 

Nl First dimension of array F 

N Number of points to be plotted 

XO First abscissa value 

DELX Abscissa increment 


Subroutine PLOTER (FX, XA, HMAX, LAMBDA, IB, NWAVE) 
The routine initializes various values required to generate printer plots and computes 
pitch-and-heave ratios. The printer plots that are generated consists of pitch-and-heave time 


histories. A description of input variables follows. 
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FX A two-dimensional array, containing time histories to be plotted 
XA Initial time 


HMAX Time-interval increment; time interval between values in FX is given by 
HMA X*PTIME 

LAMBDA Wavelength 

IB Number of values to be plotted 

NWAVE Position in FX at which wave is completely turned on 


Function RMP (T, START, RISE) 

The RMP is a function that calculates a value between O and | corresponding to time T, 
based on a straight line from time START with a value of 0 to time START plus RISE with 
a value of 1. It is used to lower the initial wave amplitude to avoid large transients at start 
of the computations. 


The arguments are: 


Wy Actual time 
START Time at which to begin the ramp from 0 to 1 
RISE Duration of rise from O to | 


The function reaches the value 1 at time START plus RISE, if the rise is 0.0, RMP will 


return a value of 0.5. 


Subroutine TRAP (F, DX, NPTS, ANS) 
This routine performs the evaluation of an integral using a trapezoidal approximation. 


The argument variables are defined as follows: 


F Array of integrand values 

DX Increments at which F is evaluated 
NPTS Number of values in F 

ANS Result, which is equal to 


NPTS 
DX bs F(i) - 0.5 [F(1) + rovers) 


i=l 


PROGRAM PLTHSP 

This program uses a data file created by program MAIN to create CALCOMP plots. The 
data are read from logical Unit 9 and are rewritten on Unit 7 for CALCOMP input. Program 
PLTHSP sets the tape output unit equal to 7 and calls SUBROUTINE CALPHI to execute the 


plot procedures. 
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Subroutine CALPLT 

This subroutine manages all the I/O operations and performs the necessary calculations 
required to generate the plots. After reading the card data (two or three cards) subroutine 
READT is called to read the data file (Tape 9) created by program MAIN. The CALCOMP 
initializing routines are called next, after which a call to subroutine ESCALE calculates the 
necessary scaling factors. Subroutine EXAXIS is called next to determine the placement of 
the plot tick marks and identifying digits. The CALCOMP plot-generation subroutines are 
now called and, depending on the option defined by the IA parameter on card 2, plots of 


pitch and heave at the bow and CG location are generated as functions of time if IA = 1. 


Subroutine EAXIS 

The subroutine is analogous to the CALCOMP AXIS routine. The only exception is that 
the tick marks are not necessarily inch, and the height of the characters is defined by the 
input parameter HT. Function NDIGIT is called to determine the number of digits necessary 


to print an even increment of the plots functions on the axis. 


Subroutine ESCALE, ADJUST, and FUNCTION UNIT 

These subroutines find the scale to be used on the plot axis. Function UNIT is called 
to determine the axis increment size after which subroutine ADJUST is called to extend the 
minimum (AMIN) and maximum (AMAX) values so that they are even multiples of the axis 


increments. 


FUNCTION NDIGIT 

This function finds the number of digits necessary to print even increments of the 
function on the axis. Both the number of places in the entire number (NDIGIT) and the 
number of decimal places (ND) are determined, after which the value of each increment on 
the axis (ANUM) is calculated. 


Subroutine READT 

This subroutine reads the data file created by program MAIN. Data file records are 
read until the message end of file is encountered. Each record is read in the same format as 
it was written in MAIN. The information is printed to allow the user to inspect the created 
file. 
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LISTING OF COMPUTER PROGRAM FOR MOTION COMPUTATIONS 


PROGRAM MAIN(INPUT »OUTPUT »s TAPES=INPUT 9 TAPE6=OUTPUT 9 TAPE3=5125 MAIN 

e TAPE2=5120TAPES=5129TAPES) MAIN 

C MAIN 
REAL IT »K »LAMBDA 9M oMA oMMAX gNoNCGgNUgMASSyNLoIAgKAR MAIN 
INTEGER END MAIN 

Cc MAIN 
DIMENSION x (6) 9FX(25400) MAIN 

Cc MAIN 
COMMON /CONST/ NCGoECGoPI sDPR»RPD »GRAVTY »RHO oK »NUM9MA(120) »CDoTAy MAIN 

e B(120) sBETAsHW (120) »TZsVRAGoWoXDoT oXPoMolITs MAIN 

a DELTAS o TXoEST (120) »CoROsKARgMMAX (1 0) 9 TEST (120) 9 MAIN 

6 N(120) »PHALF MAIN 

* COMMON ASHIP/ MASS oCINT sQA9CE 9CE29CE3oDMU sEDMU »E2DMU 9 E3DMU 9 BF 9 BMM g MAIN 

6 NL oFLoIAoE (120) MAIN 
COMMON /IN/ BM(120) 581 (120) »VELIN MAIN 
COMMON/UUT /NPRINT »NPLOT 9 END MAIN 
COMMON/TERMS/T1 o T20T39T4oTS 9 T6oT7 9718 MAIN 
COMMON /SEAWAVE/ START »RISE,RAMP MAIN 
COMMON /INTER/ II oKTT(10) sDIFF (10) MAIN 
COMMON /IN2/ NO(120) »XAoXEoHMAX gHMIN9 A (6) 9EPSE (6) »LAMBDA MAIN 
COMMON /ACCEL / XACCL»BWACL »CGACL »BL MAIN 

Cc MAIN 
CALL INPUT MAIN 

Cc MAIN 
Cc COMPUTE INTEGRATION INTERVAL INFORMATION MAIN 
( MAIN 
NLESS = NUM-1] MAIN 

I=l MAIN 

II = 1 MAIN 
DIFFER = EST(1el) -EST(I) MAIN 
KTT(II) = 1 MAIN 

OIFF (II) = DIFFER MAIN 

DO 25 I=2eNLESS MAIN 
DIFFER= EST(1+1)-EST(I) MAIN 
KTT(II) = KTT(II) 1 MAIN 

IF (DIFFER NE UIFF(II))GO TO 24 MAIN 

GO TO 25 MAIN 

24 II = II*l MAIN 
KTT(II) = 1 MAIN 

DIFF (II) = DIFFER MAIN 

25 CONTINUE MAIN 
KTT(II) = KTT(II) ol MAIN 

Cc # * # # # CHECK IF NUMBER OF INTERVALS EXCEEDS DIMENSION MAIN 
IF (11eGT.10) WRITE (6928) (KTT(I) sDIFF (I) sI=loII) MAIN 

IF (IIeGT.10) STOP 4 MAIN 

C # * # # # POINT AT WHICH MULTIPLE RUNS START MAIN 
8 CONTINUE MAIN 
TIME=XA MAIN 
KOUNT=1 MAIN 
END=END=-1 MAIN 

WRITE (6939) MAIN 

39 FORMAT (1H1) MAIN 
c# *# # # # # # # READ IN INITIAL CONDITIONS MAIN 
Cc X(1) = VELOCITY» X(2) = 2 DOTs X(3) = THETA DOT MAIN 
Cc KX(4) = Xo X(5) = Z9 X(6) = THETA MAIN 
Cc THETA IS pPEAD IN DEGREES THEN CONVERTED TO RADIANS IN PROGRAM MAIN 
C MAIN 
READ (5510) (X(1) 9 I=1 56) MAIN 

Cc MAIN 
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Cc DATA » USED IN RAMP FUNCTION, TO TURN ON WAVE MAIN 
READ (5510)START sRISE MAIN 

Cc MAIN 
30 FORMAT (8F 19.4) MAIN 
C# * # © # @ # » WRITE OUT THE INPUT VALUES MAIN 
WRITE (6019) START»RISE KAR MAIN 

19 FORMAT (" START = "FIO 04% 9/5" RISE = "oF100% 9/9" KAR = "oFLOMAIN 
04) MAIN 

Cc MAIN 
Cc TME IS THE TIME AT WHICH THE INTEGRATION INTERVAL IS MAIN 
c TO KE CHANGED MAIN 
c HMX IS THE NEW MAXIMUM INTERVAL SIZE AFTER TIME TME MAIN 
c HMN IS THE NEW MINIMUM INTERVAL SIZE FOR KUTMER TO SUB-DIVIDE MAIN 
Cc THE MAXIMUM INTERVAL UP TO MAIN 
Cc IF THIS OPTION IS NOT USED SET TME TO THE STOP TIME OF THE RUN MAIN 
Cc MAIN 
READ (5510) TME »HMX»HMN MAIN 

WRITE (6911) TME »HMAX »HMX 9 HMINoHMN MAIN 

ll FORMAT(® AT TIME #9F7.29% THE MAXIMUM INTERVAL SIZE FOR INTEGRATIMAIN 
#ON WILL BE CHANGED FROM oF 10049% TO @9Fl0.49/5 MAIN 

#@ AND THE MINIMUM SIZE FOR HALVING CHANGES FROM #5F10,4, MAIN 

# ® TO oF 10.4) MAIN 

C ADJUST THE TIME FOR CHANGE OF INTEGRATION INTERVAL MAIN 
Cc FOR CHECK AGAINST TIME IN THE INTEGARTION LOOP MAIN 
T™ = TME-(HMAX/2.) MAIN 

Cc SET SWITCH FOR CALCULATION UF PITCH AND HEAVE RATIOES MAIN 
Cc ON NEXT CALL TO PLOTER MAIN 
IPT = 0 MAIN 

IF (TME.EQ.XE) IPT = 1 MAIN 

Cc MAIN 
READ(5,10) PERCNT MAIN 

XACCL = ECG=PERCNT#BL MAIN 

WRITE (6912) PERCNT»XACCL MAIN 

12 FORMAT(*® THE x USED FOR THE BOw AND CG ACCELERATION COMPUTATIONS MAIN 
oIS EQUAL TO ECG-*sF10.497H*BL OR 0F10.4) MAIN 

C MAIN 
WRITE (6023) MAIN 

WRITE (6947) MAIN 

23 FORMAT(1H ,//) MAIN 

47 FORMAT (" STATION NO. 53Xo"DEAD RISE" 9 BX o"EST"58X 9 "NOM"s MAIN 

# 1OX»"BEAM") MAIN 
WRITE (6955) ((I1sBETASEST (I) »NO(T) 9BM( 1) ) 9 T=) 9 NUM) MAIN 

55 FORMAT (OXoT20SX oF 10049 4XoF1004%94X9F 10040 3KX0F 1004) MAIN 
WRITE (6923) MAIN 

WRITE (6956) (X(I) oI=196) MAIN 

56 FORMAT(" X VALUES" 54X96(F10,492X)) MAIN 

Cc # # # # # @ @ # CHANGE INPUT FROM DEGREES TO RADIANS MAIN 
X(3) = X(3)#RPD MAIN 

X(6) = X(6)#RPD MAIN 

c MAIN 
WAVE = STA2T+¢RISE MAIN 

NWAVE = 0 MAIN 

C# # # # #@ # # & WRITE OUT COMPUTED ARRAYS MAIN 
WRITE (6957)Mo IT 9K oC oPHALF oP I »GRAVTY MAIN 

IF (NPRINT.LT.%) GO TO 62 MAIN 

WRITE (6958) (E (1) oT=1yNUM) MAIN 

WRITE (6559) (N(I) oI=1 9NUM) MAIN 

WRITE (6064) (MMAX(I) » I=] 9NUM) MAIN 

WRITE (6065) (TEST(I) 9 I=) 9NUM) MAIN 
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101 
102 
103 
104 
105 
106 
107 
108 
109 
110 
111 
112 
113 
114 
11S 
116 
117 
118 
119 


62 
28 


CONTINUE 
WRITE (6928) (KTT(I) oOIFF (I) sT=)e II) 
FORMAT(#® KTTsDIFF #511092X,F10,4) 


MAIN 
MAIN 
MAIN 


ST FORMAT (4H M= 9F10.494H T= oF 10.4 94H K= 9F100494H C= 9F10.4,11H PI#MAIN 


Gre® 


aaN0 


asl 
c*#* 


Go 


98 


c## 
99 


@RHO/2= 9F10.495H PI= 5F10.4,10H GRAVITY= 5F10.4) 


FORMAT (" F(1)%510F1004) 
FORMAT (" N(I)% 9 10F10.4) 
FORMAT ('* MMAX(1) "5 lOF 1004) 
FORMAT ( TEST(1)%sl0F 10.4%) 
IB = 1 

IPRINT = NPRINT 

WRITE (4991) 


# @ # @ # # WRITE HEADINGS AND CONDITIONS AT TIME = 0. 
Ql FORMAT (LH1 2X o"T IME s9X9"XDOT" o9Xo"ZDOT" s9Xo"THETA DOT's 6X9 


6 LAX99Xo1HZ99X o9SHTHE TA 9X a QHNL 9 9X9 2HFL 
@ 4Mo8HBOW ACCL>4X97HCG ACCLo//) 

WRITE (4992) TIME» (X(I) sIT=1 96) sNL FL »BWACL »CGACL 
WRITE(9) TIME» (X(1) 5 1=496) »BWACL »CGACL 
KOUNT = KOUNT?1 

FX(101B)=X(5) 

FX (291B) =X (6) 

IKUTM= (XE=xA) /HMAX*.05 

IKUTM = (TME=XA) ZHMAX @ (AE=TME) SHMX @ .05 
FIRST=0.0 

NEQS=6 

IKUTS=0 


START OF INTEGRATION LOUP 


CONTINUE 
NPRINT = IPRINT 
# # @ # # ® CHECK PITCH .GT, 5236 RADIANS 
IF (X(6) eGT.05236)G0 TO 853 
® ® @ # © » PERFORM INTEGRATIONS 
IF (TIME LT. TMeOReTMEEQeXE) GU TO 98 
IF (IPT.EQ.1) GO TO 98 


HMIN = HMN 

HMAX = HMX 

FIRST = 0.0 
CONTINUE 


CALL KUTME9 (NEQS 9 TIME »HMAX 9X 9EPSE sAgHMINgF IRST) 
IKUTS=IKUTS¢1 

IF (FIRST.E9.2)GO0 TO 861 

IF (KOUNT NF el eANDeKOUNTeNE 041) GO TO 99 

WRITE (4991) 

KOUNT=1 

*# # @ # # # WRITE OUT TIME INTERVAL RESULTS 
WRITE (4992) TIME» (X(I) 9 I=] 96) s9NLoF lL oBWACL »CGACL 
WRITE (6993) Tl oT2oT39T4eTSoT65T7 9 T898MMoBF 
WRITE(9) TIMEs (X(1) 9 1=496) »BWACL »CGACL 

IF (TIME sLT.TMceOReoTMEcEQeXE) GU TO 200 

IF (IPT.EQe1) GO TO 200 

CALL PLUTE? (FX5XAgHMAX »LAMBDA sIBoNWAVEs IPT) 

IPT = 1 
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175 
176 
177 
178 


NAAANANANANANAANAANAANAANAINAND 


20 CONTINUE 


1821861 


FX(1,1B)=x(5) 
FX(2s1B) =x (6) 
93 FORMAT(" ",10E10.4) 
92 FORMAT (1Xo11(F10.402Xx)) 
190 CONTINUE 
~ KOUNT=SKOUNT ¢1 


IF (NWAVE.GT.0)GO TO 21 
IF (TIME GT. WAVE) NWAVE=KOUNT 


21 CONTINUE 


IF (TIME LE. XE AND. IKUTSoLT.IKUTM)GU TO 851 


WRITE (29852) 

854 CONTINUE 

852 FORMAT( “ END OF KUTMER") 

853 CONTINUE 
CALL PLUTE2 (FX »XAoHMAX »LAMBDA s IB yNWAVE 5 IPT) 

# # @ @ # @ # CHECK FOR LAST RUN IF NOT CYCLE HACK TO READ 
NEW DATA FOR NEXT RUN 


IF (ENO.NE.1)GU0 TO 8 


GO TO 999 
# # @ # KUTMER ERROR MESSAGES 


861 WRITE (69862) 
862 FORMAT ("' ERROR CRITERION IN KUTMER CAN NOT BE MET") 
WRITE (6956) (X(I) 9I=196) 
WRITE (6986) TIME 
86 FORMAT (" TIME ='"9F10.4) 
IF (END.NE.1)GO TO 8 
GO TO 853 
999 CONTINUE 
END FILE 9 
END 
SUBROUTINE PLUT2 (F oF MINoF MAX gNVAR »NFUN NI 9N9 XO 9DELX) 
PLUT FIRST N POINTS OF UP TO 26 FUNCTIONS F(X) 
F(I9J) CONTAINS THE VALUE FOR THE JTH POINT OF THE ITH FUNCTION 
FMIN(I) AND FMAX(I) CONTAIN THE MIN AND MAX ORDINATE VALUES FOR 
THE ITH FUNCTION, 
NVAR(I) AN ARRAY OF TITLES FUR THE VARIOUS FUNCTIONS 
TO BE PLOTTED AGAINST THE ABSCISSA 
NF UN NUMBER OF FUNCTIONS TU BE PLOTTED - OIMENSION OF 
NVAR» FMIN» FMAX 
Nl USED ONLY IN F(Nl»l) AS PASSED DIMENSION 
N NUMBER OF POINTS IN A SINGLE PLOT FRAME 
x0 FIRST ABSCISSA VALUE 
DELX ABSCISSA INCREMENT 


1 


2 
ro 


DIMENSION 2STEP (26) oF (NI oN) oF MIN (NFUN) oF MAX (NFUN) » VLAST (26) » 
VF 19ST (26) HEAD (6) »STEP (26) 
INTEGER CH(26) »NVAR( NFUN) »DOT»ASTERsPLUS » BLANK 
INTEGER C 
INTEGER A (101) 


DATA BLANK ,DOT,ASTER»PLUS/IH olHeolH®slH+e/ 

DATA CH(1) »CH(2) »CH(3) oCH (4) »CH(S) »CH(6) »CH(7) »9CH(8) »CH(9) »CH(10) 
7 HA 4» LHB » LHC » LHD » LHE » LHF » 14G » IHH oIHI olHJ # 

DATA CH(11) »CH(12) »CH(13) 9CH(14) »CH(15) »CH(16) »CH(17) sCH(18) 

f 7 AK 4 LAL » JAM » JAN » JHO » LHP » ILH@ » ILHRS 

DATA CH(19) »CH(20) oCH(21) »CH (22) »CH(23) »CH (24) »CH(25) »CH (26) 
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2 


4 HS » HT » IHU » LHV 5» LHW 9» 


IF (NFUN.LE.0.OR.N.LE.0) RETURN 


C PRINT HEADINGS, 


46 


30 


WRITE (6946) 


FORMAT (///) 


DO 40 I=l 5NFUN 


TENM=A0S (FMAX(1) <FMIN(I)) 
EXP=1l. 
IF (TENM.EQ.0-) GO TO 2 


C ARING TENM TO A VALUE BETWEEN 1 AND 10 


3 


C SET 


IF (TENM.LTo1le) GO TO 1 
IF (TENM.LT.19.) GO TO 2 
EXP=EXF#1", 
TENM=TENM*# .1 

GO TO 3 

EXP=EXP#,1 
TENM=TENM#10.6 

IF (TENM.GT.1.) GO TO 2 
GO TO 1 

UP VALUE SETWEEN GRID LINES, RSTEP. 
PSTEP=5, 

IF (TENM.GE.5.)PSTEP=10,. 
IF (TENM.LT220e)PSTEP=2.6 
RSTEP (1) =PSTEP#EXP#,} 


C CUMPUTE VALUE OF STARTING LINE» VFIRST.o 


FIRST=FMIN(I) /RSTEP(I) 

IF (FMIN(I) LT O-)FIRST=FIRST=-1l, 
FIRST=AINT(FIRST) 

VF IRST(1)=FIRST#RSTEP (I) 


C CHECK END LINE VALUE» VLAST. 


VLAST(I)=VFIRST (1) +10.#RSTEP (1) 
IF (VLAST (I) oGTeFMAX(I))GO TO & 


C IF GRAPH IS TUO SMALL TAKE NEXT LARGER STEP. 


AA=PSTEP 

IF (AA.LT.5.)PSTEP=5, 
IF (AAcEQe5e) PSTEP=106 
IF (AA.LT.2100) GO TO S 
PSTEP=e. 

EXP=10.#ExXP 

GO 70 5 


C CUMPUTE VALUE BETWEEN POINTS»STEP. 


4 


6 
40 


45 FORMAT (1X9A193H = 9Al095X91PE12.495(8X91PE1204) ) 
DO SO J=19101 


50 
55 


STEP (I) =RSTEP (1) #1 
RK=0. 


DO 6 KK=196 


HEAD (KK) =VFIRST(I) +2. #RK*®RSTEP (1) 
RK=RKe1,. 


LHX 95 


WRITE (6945) CH(I1)» NVAR(I)»s (HEAD(KK) oKK=1 96) 


A(J) =BLANK 

IF (MOD (Jo1%) eEQel) A(J)200T 

CONTINUE 

WRITE (6555) AoA 

FORMAT (25Xo1l01A1/15X 54HTIME s6X59101A1) 


C PLUT EACH PUINT 


DO 100 J=1,N 
B=X0+FLUAT (J=-1) *DELX 
DO 70 K=19101 
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70 


4@0 


39 


A(K) =BLANK 
IF (MOD (K910) EGol) A(K) 200T 
IF (MOD (J9S )cEQel) A(K)=D0T 
CONTINUE 
DO 80 I=1,NFUN 
LOC=((F (1, J) -VFIRST(I)) STEP (1) 1.5) 
C=A(LOC) 
A(LOC) =CH(T) 
IF (C.NEeBLANKeAND.C.NE.DOT) A(LOC)=ASTER 
CONTINUE 
IF (MOD (J910).EQ.1)GO TO 95 
WRITE (6985) A 
FORMAT (25X9101A1) 
GO TO 100 
WRITE (6915)BoA 
FORMAT (12X1PEL2.451Xo101A1) 
CONTINUE 
RETURN 
END 


SUBROUTINE KUTMER(ND9T oHo V0 ,EPSE sAgHCX oF IRST) 
DIMENSION Y0(6) sV1 (6) 9¥2(6) 9F0(6) 9F 1 (6) oF 2 (6) g9EPSE (6) 9A (6) 


COMMON /O0UT /NPRINT oNPLOT END 
COMMON /ACCEL / XACCLoBWACL »CGACL oBL 
DATA NAM1L»NAM2 /2HY1ls2HY2 / 


ND = NUMBER OF EQUATIONS, NO, UF COMPONENTS OF YO 

T = INDEPENDENT VARIABLE 

H = INCREMENT FOR WHICH SULUTION IS TO BE RETURNED + OR = 
yO = THE VECTOR OF DEPENDENT VARIABLES. ENTER WITH INITIAL 


VALUES AT T AND RETURN WITH VALUES AT TeH 


EPSE = RELATIVE ERROR CRITERION FUR COMPONENTS OF YO .GT ABS(A) 
A = ABSOLUTE ERROR SRITERION FUR COMPONENTS OF YO .LT. ABS(A) 
NUTE-- EPSE AND A MUST BE SPECIFIED FOR EACH COMPONENT OF THE SYSTEM 


HCX = THE SMALLEST STEP SIZE USED IN THE INTEGRATION 


IF Iv IS ENTERED WITH A CHANGED H 


IF (FIRST) 20910520 
eee e ese oe FIRST ENTRY 


ees 227 S20 OTHER ENTRY. 


r 
o 
(o) 
(—) 


HCX HC 
IF (HC.NE.9.) GU TO 30 
WRITE (69800) 


FIRST = 2. 
RETURN 
e eee e ce = S CALLS TO DAUX 


0 CALL DAUX(T»Y09F0) 


IF (NPRINT.EQo5) WRITE (69400) YOoT»FO 
FORMAT (6 (2X 9F 10.4%) o4HTIME 92X oF 10.4%) 
IF (NPRINT.F£Q.5) WRITE (69400)HC 

DO 40 I=1sND 


D5) 


FIRST SHOULD BE 0 WHEN KUTMER 1S ENTERED FOR THE FIRST TIME 
AFTER THAT FIRST IS 1 IF KUTMER IS ENTERED WITH THE SAME H OR 


0 FORMAT (SXo4S5HKUTMER ENTERED WITH ZERO INTEGRATION INTERVAL ) 


PLOT2 88 
PLOT2 89 
PLoT2 90 
PLOTe 91 
PLOT2 92 
PLOT2 93 
PLOT2 94 
PLOT2 95 
PLOT2 96 
PLOT2 97 
PLOT2 98 
PLOT2 99 
PLOT2100 
PLOT2101 
PLOT2102 
PLOT2103 
PLOT2104 
PLOT2105 
PLOT2106 
KUTMER 
KUTMER 
KUTMER 
KUTMER 
KUTMER 
KUTMER 
KUTMER 
KUTMER 
KUTMERI1O 
KUTMERI1 
KUTMERI2 
KUTMERI13 
KUTMERIS 
KUTMERIS 
KUTMERI6 
KUTMERI7 
KUTMERI8 
KUTMERIO 


ODBNANFSF WN 


IF FIRST IS 2 THE ERROR SRITERIA CANNOT BE MEET AND THE STEP SIZE IKUTMER20 
REDUCED TO 47128. 


KUTMER21 
KUTMER22 
KUTMER23 
KUTMER24 
KUTMER2S 
KUTMER26 
KUTMER27 
KUTMER28 
KUTMER29 
KUTMER30 
KUTMER31 
KUTMER32 
KUTMER33 
KUTMER34 
KUTMER3S 
KUTMER36 
KUTMER37 
KUTMER38 
KUTMER39 
KUTMER4S0 
KUTMER41 


40 
Cc 

50 
Cc 

60 
Cc 

70 
Cc 

80 
Gien= 
Cee 

85 
Giles 

87 
Creu 
Cee 

90 
Cumin 

91 

92 
Cee 

95 


710 


720 


VICI) = YO (1) *(HC/3.) #FO(1) KUTMERS2 
IF (NPRINT .€Q.5) WRITE (65400) YI oT KUTMER43 
KUTMER44 

CALL DAUX(T*HC/3¢ V1 oF 1) KUTMER4S 
IF (NPRINT o£ Qo5) WRITE (69400)F 1 9T KUTMER4G6 
Q0 50 I=l»ND . KUTMER47 
YI(I) = YO(1) *(HC/6.) #FO(1) + (HC/60) #F 1 (1) KUTMER48 
IF (NPRINT 0 Q05) WRITE (65400) V1 oT KUTMERS9 
KUTMERSO 

CALL DAUX(T*¢HC/3.9YV19F 1) KUTMERS1 
IF (NPRINT o=Q05) WRITE (69400)F19T KUTMERS2 
DO 60 I=l»ND KUTMERS3 
YI(I) = YO(1)+(HC/8.) #F0(1) +. 375#HC#F] (I) KUTMERSG 
IF (NPRINT £EQ05) WRITE (65400) VI oT KUTMERSS 
KUTMERS6 

CALL DAUX(T*HC/269VlsF2) KUTMERS7 
IF (NPRINT o=Q05) WRITE (69400)F20T KUTMERS8 
00 70 I=l»ND KUTMERSS 
YL(T) = vO(1)*(HC/2.) #F0(1) =) SeHCeFl (1) +2, ehnC#Fe (1) KUTMER60 
IF (NPRINT 5=Q05) WRITE (65400) V1 5T KUTMER61 
KUTMER62 

CALL DAUX(T+HCoV1 oF 1) KUTMER63 
IF (NPRINT .=Q05) WRITE (69400)F 1 9T KUTMER64 
DO 80 I=1l>5ND KUTMER65 
Y2(1) = YO(I) *HC/6.#F0(1) + (20/36) *HC#F2 (1) + (HC/6.) #F 1 (T) KUTMER66 
IF (NPRINT .£Q.5) WRITE (69400) Y20T KUTMER67 
INC = 0 KUTMER68 
ese eer eee CHECK ERROR CRITERIA KUTMER69 
DO 110 1=1.NO KUTMER7O 
ZZZ = ABS(YIL(I))=-A(I) KUTMERT1 
IF (222) 85587,87 KUTMERT2 
soso = 2 = = ABSOLUTE ERRUR KUTMER73 
ERROR = ABS(.2#(Y1(I)-Y2(1))) KUTMERT4 
IF (ERROR=A(1)) 100,100,590 KUTMERTS 
cee cee eee RELATIVE ERKOR KUTMER76 
ERROR = ABS(e2-e2"Y2(1)/VI(1)) KUTMERT7 
IF (ERROR=EPSE(I)) 1009100990 KUTMER78 
Fe a ie Sule SINCE ERROR eGT. ERROR CRITERIA CHECK IF HC.GT H/KUTMERT9 
ee eee ee IF YES THEN HALVE INTERVAL. OTHERWISE STOP, KUTMERBO 
X = 128.#AS (HC) =ABS (H) KUTMERB1 
IF (X) 91,95,95 KUTMERB2 
-2f2ee-+-ee ERROR TOU LARGE KUTMER83 
WRITE (6992) I9T sERROR®HC KUTMER84 
FORMAT (/18H FOR EQUATION NO. I2927Hs THE RELATIVE ERROR AT T = KUTMERBS 
e E15e8>5 4H IS 9£15.8913H STEP SIZE = 5615-8) KUTMER86 
FIRST = 2. KUTMERB7 
RETURN KUTMER88 
=<s-e--- = HALVE INTERVAL KUTMERB9 
HC = HC/2, KUTMER9O 
IPLOC = 2#IPLUC KUTMERS1 
LOC = 2*LUC KUTMER92 
HCX = HC KUTMER93 
WRITE (29719) Ts I s9ERROR9HC KUTMER94 
FORMAT (/8H TIME = 9Ff10¢395Xs26HHALVE INTERVAL. EQUATION 9135 KUTMERSS 
6ol3H HAS ERDOR = 5£16.896X917H STEP SIZE NOW = 5615.8) KUTMER96 
WRITE (29729) NAM2 9 (Y¥2(J) 9 J=1 ND) KUTMEROT 
WRITE (29720) NAM» (Y1 (JU) 9 J=19ND) KUTMER98 
0 FORMAT( 2XeA2 / 3(10E1305/)) KUTMER99 
Go TO 390 KUTME100 


56 


Cer eee22-222 2 TEST IF INTERVAL LENGTH CAN BE DOUBLED 


KUTME1O] 


100 IF (ERRUR#64.-EPSE(I)) 11051205101 KUTME102 
lol INC = 1 KUTME103 
110 CONTINUE KUTME 104 
C-- 2-2 - == - = = UPDATE T AND SOLUTION KUTME105 
ll T = TeHC KUTME106 
DO 112 I=1,ND KUTME107 

le YO(I) = v2(1) KUTME 108 
Ce - ee ee = = - = ~ GET SOLUTION IN NEXT INTERVAL KUTME109 
LOC = LUC+1 KUTME110 

IF (LOC-IPLOC) 12092105210 KUTME LLL 

120 IF (INC)210,130,210 KUTME112 
130 IF (LOC-(LOC/2) #2) 21051405210 KUTME113 
146 IF (IPLOC-1) 210,210,200 KUTME114 
C--+22e-+2+- - = = DOUBLE INTERVAL LENGTH KUTME115 
200 HC = 2,%HC KUTME116 
LOC = LUC 72 KUTME117 

IPLOC = IPLOC/2 KUTME118 

210 IF(IPLOC-LOC) 30,329,30 KUTME119 
329 BwACL = F0(2)-xACCL#F0(3) KUTME120 
CGACL = FO0(2) KUTME121 
RETURN KUTME 122 

END KUTME123 

END KUTME124 
SUBROUTINE DAUX (TIME »X»RHS) DAUX 2 

Cc DAUX 3 
C TIME TIME AT WHICH SYSTEM IS TO BE EVALUATED DAUX 4 
© x STATE VECTOR DAUX 5 
Cc RHS THE RIGHT HAND SIDE UF THE EQUATION S =F A DAUX 6 
Cc DAUX 7 
REAL KAR DAUX 8 

REAL IAs IT oMoK sMAgMASSoNCGoNL oN gMMAX DAUX 9 
INTEGER END »PTIME DAUX 10 
DIMENSIUN x (6) »RHS(6) oF (391) 9A(393) 9 INDEX (353) 9 DAUX 11 

Aj R(120) 0V(120) 9D(120) DAUX 12 

Cc DAUX 13 
COMMON /SHIP/ MASS oCINT »QA9CE 9CE2 9CE3 yDMUsEDMU,E2DMUsE30MU 9 BF 9BMMyDAUX 14 
NLoFLsIAgE (120) DAUX 15 

“COMMON /CUNST/ NCGsECG oP sDPR»RPD sGRAVTY »RHOsK »NUM»MA(120) »CDy9 TAs OAUX 16 

e B(120) sBETA HW (120) 9 TZ 9DRAGoWoXDoT oXPoMolTs DAUX 17 

0 DELTAS s TX9EST (120) »CoRO»KAR gMMAX (1-0) oe TEST (120) o DAUX 18 

N (120) »>PHALF DAUX 19 

COMMON /IN/ BM (120) ,81 (120) »VELIN DAUX 20 
COMMON/UUT /NPRINT »NPLOT sEND VAUX 21 
COMMON /SEAWAVE/ START »RISE ,RAMP DAUX 22 
COMMON /WAVE/ RoPT(120) oZMAgZWMA EMAS 5 ZZWMA 9 ZWEMA 9 ZOWMA 9ECMAZ 9 DAUX 23 

6 ZwOOT (120) DAUX 24 

Cc DAUX 25 
RAMP = RMP(TIME START RISE) DAUX 26 

PIH = PI/2, DAUX 27 

CT = C#TIME DAUX 28 

Cx6 = CUS(x(6)) DAUX 29 

SK6 = SIN(X(6)) DAUX 30 
CaeeeeeeSET VALUFS UF MA AND B DAUX 31 
00 75 I=1,NUM DAUX 32 

PT(L) = (X(4) E (1) #CX6eN(1) #SA66CT) HK DAUX 33 

R(I) = RU®COS(PT(I)) #RAMP DAUX 34 

Cc # * # # # # # # COMPUTE HW SUBMERGENCE OF A POINT AND R- THE WAVE DAUX 35 
Cc HW(T) IS IN THE FIXED COORUINATE SYSTEM DAUX 36 


Si 


Cc 
65 

Cc 

Cc 

Cc 

Cc 

Cc 
70 
75 
74 
76 
77 
78 
79 
80 
si 
82 
85 

Cc 

Cc # # 

Cc 
15 
17 

Cc # * 
18 

C # # 


HW(I) = X(5)-E (1) #SX6eN(1) #CX6-R (1) 
IF (HW(1) GT.) GO TO 65 
CRAFT IS NOT SUBMERGED 
MA(I) = 0. 
Bl(1I)=0. 
B(I) = 0. 
GO TO 75 
V(I) -RO#K*SIN(PT (1) ) *RAMP 
0 (1) Hw(I)7(CA6—-V (1) #SX6) 
D(I) 1S IN THE BODY AAIS SYSTEM AND IS THE SUBMERGENCE 
IF (D(1) GE.TEST(I)) GO TU 70 
CRAFT IS PARTLY SUBMERGED 
B(I) = O(1)#(1./TA)#PIH 
BL(T> D(T)#(1-/TA) #PI 
MA(I) KAR#PHALF #B (I) #B(1) 
GO TO 75 


CHINE IS IMMERSED 
Bl ARRAY IS USED FOR THE INTEGRALS OVER THE PORTION 
OF THE HULL FOR WHICH THE CHINE IS NOT IMMERSED 

MA(I)=MMAX(1) 
B(I)=BM(1) 
B1l(I)=06 
CONT INU 
IF (NPRINT.LTe%) GO TO 85 
WRITE (6974) TIME 
FORMAT(" TIME = "yF10.4) 
WRITE (6976) (X(T) » T=1 96) 
WRITE (6977) (R(I) » I=] NUM) 
WRITE (6978) (HW(I) »I=1 NUM) 
WRITE (6079) ( B(I) »I=l NUM) 
WRITE (658%) (V(1I) o I= »9NUM) 
WRITE (6581) (O(1) » l=] NUM) 
WRITE (6982) (MA(I) » T=) »NUM) 
FORMAT(" X(I) "“96(2X0E120c6)) 
FORMAT (* R(T) %910F 10.4) 
FORMAT ( HW(I)%510F10.4) 
FORMAT (" 8(1)%510F10.4) 
FORMAT (* v(I)%510F10.4) 
FORMAT (" (I) %910F 10.4) 
FORMAT (" MA(I) “slOF1064) 
CONTINUE 


* @ # # # » COMPUTES NL AND FL ANU THE ASSOCIATED INTERGALS 
CALL FUNCT (X) 


IF (NPRINT.LT2%)GO TO 17 

WRITE (6915) TXoFLoDRAGoTZoWgNboXDoT oXP 
FORMAT(" ",10E12.6) 

CONTINUE 

# # # # # COMPUTE THE F VECTOR 
F (lol) = \TX#FL#SX6=DRAG#CX6 
F(151)=0.0 

F(20l) = T2*FL#CX6*DRAG*SX6eW 
F (391) =NL=DRAG#XD+T#xXP 

IF (NPRINT.LT.3)G0 TO 18 

WRITE (6910) (F (Tol) »I=193) 
CONTINUE 

# &# # # # COMPUTE THE A MATRIX 
A(lol) = M+MASS#SX6#SX6 


58 


DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUXK 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
VAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
OAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 


A(lo2) = MASS*®SKX6#CX6 
A(leo3) = =)A*#SX6 
A(ls2) = oe 

A(l9o3) = 


Al2el)= AIRE) 

A(202) = M+MASS#CX6#CX6 

A(293) = -QA#CX6 

A(391)=A(1,.3) 

A(392)=A(2.3) 

A(393)=SIT*IA 

IF (NPRINT.LT2£3)GO TO 25 

WRITE (6012) (A(I 91) oI=193) 

WRITE (6913) (A(I 92) oL=1 93) 

WRITE (6914) (A(1,3) »I=1,3) 

# @ # @ # INVERT THE A MATRIX 

CALL MATINS (A93939F 91919DETERMsIDo INDEX) 
IF (ID.EQ.2) WRITE (6026) 

FORMAT (*! MATRIX IS SINGULAR my) 


Ceeee#eA ON RETURN WILL CONTAIN THE INVERSE MATRIX 


c*#* 
25 
26 

Cc 

Cc 

Cc 

c## 
10 
12 
13 
14% 
39 
35 
40 

Cc 


I0=2 MATRIX IS SINGULAR 
=1 INVERSE WAS FOUND 


® # # © # COMPUTE THE RIGHT HAND SIDE 


RHS(1) = F(lsl) 
RHS(2) = F(291) 
RHS(3) = F(391) 
RHS(1) = 0.0 
RHS(4) = X(1) 
RHS(5) = X(2) 
RHS(6) = X(3) 


FORMAT (" F(T oh) %93(2X5E12.4)) 
FORMAT ("' A(Tol) "s3(2XsE12.4)) 
FORMAT (" A(T92) "s3(2XsE12.4%)) 
FORMAT ("" A(1593) "93(2XsE1264)) 
IF (NPRINT.LT.2) GO TO 40 

WRITE (6912) (A(Io1) ,IT=193) 

WRITE (6013) (A(1,92) »I=1 93) 

WRITE (6914) (A(I93) 9I=193) 

WRITE (6935) (RHS(I1) »I=1 96) 
FORMAT ("* RHS(I) %96(2X9E1206)) 
CONTINUE 

RETURN 

END 

SUBROUTINE FUNCT (X) 

REAL KAR 


REAL IAs TAAgIPART 9K oKPI o>MAgMASS NL »NCGo IT oMoMMAX oN 


INTEGER END 
DIMENSION IPART (120) C1 (120) »C2(120) 5 


26 (120) oZ7(120) 
0X (6) 9 VMAA (120) 


e@ee@fFf 


11 (120) 902 (120) 903 (120) 904(120) 505 (120) 906(120) » 
QPART (120) 9Z1 (120) 922 (126) 923 (120) 924 (120) 925 (120) » 


DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
DAUX 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 


COMMON /SHIP/ MASSsCINT sQA9CEsCE22CE3sDMU sEDMUsE2DMU 9 E 3DMUs BF s BMM oF UNCT 


NL oFLoTAsE (120) 


* COMMON JCONST/ NCGoECGoPI »DPR»oRPD »sGRAVTY »RHO 9K »NUM9MA (120) »CDo TAs 
e B(120) sBETAsHW (120) 9» TZ9DRAGgWoXDoTs XPoMolTo 
DELTAS o TXoEST (120) sCoROsKARgMMAX(1 0) op TEST (120) » 


e N(120) »>PHALF 
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FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 


& 


50 
si 


60 


COMMON /INZ BM(120) oB1 (120) »VELIN 
COMMON /U0UT /NPRINT »NPLOT »END 


FUNCT 
FUNCT 


COMMON /WAVE/ R(120) sPT (120) o> 2MAgZWMAgEMAS,ZZWMA, ZWEMA, Z2WMA,E2MAZF UNCT 


9ZWOOT (120) 


e 
COMMON /INTER/ II oKTT (10) oDIFF (10) 


COMMON /SEAWAVE/ START »RISE,RAMP 
COMMON /TEST/ VMA 
#2 2 #@ #@ » # @ INITIALIZE INTEGRAL SUMS 


0.0 

(1) ®SIN (XK (6)) oX (2) @CUS (X (6) ) 
Sx6 = SIN(x(6)) 

CX6 = COS(X(6)) 

WO = K#C 

# @ # @ # # SET UP THE FUNCTIUNS FOR THE INTEGRALS 
DO 90 I=1l»NUM 

IPART (1) =E (1) #E (1) #MA(T) 

QPART (I) =E (1) #MA(T) 

ZWOOT(I) = -RU#WO#SIN(PT (I) ) #RAMP 

U = X (1) #CX6—-X (2) *5X6e¢ZWDOT (1) #SX6 
VEL = VPART=X(3) #E (1) -ZwDOT (1) #CX6 


N 
N 
a 
= 
> 
1 0 tb 


Z1(1T) = MA(I)*#ZWOOT(T) 

Z2(1) = -Ma(1I)*#COS(PT(I)) #RAMP 
Z3(1) = E(1) #221) 

24(1) = E(T) #2101) 

Z5(1) = U#z2(1) 

Z6(1) = E(7)#25(1) 

Z7(1) = MA(I) #VEL#U 


IF (VELeLE.0.) GO TO 60 
IF (BI(1) LE. 0.0) GO TO 50 
DRDT = ZwOOT (1) #(X(L) Cox (3) ®(N( 1) @CX6=E (1) #SX6) ) /C 


DI(1) = — VEL®BL (1) #(X (2) x (3) # (CXO*E (1) #SX6#N(T)) 
GO TO 51 

DI(1) = 0. 

CONTINUE 

p2¢1) = €(1)*D1(1) 
Cl(1) = VEL#VEL#B(I) 
C2(I) = €(1)#C1 (1) 
GO TO 6 

OAS) en Oa 

D2(I) = 0. 

C1(1) = 06 

C2(1) = 0. 
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FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
(PAGE 4 OF NOFUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 


=-DRDT) FUNCT 


FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 
FUNCT 


C 


Cc 


6l 


CONTINUE 


03(1) 
04 (1) 
PIH = 


-OS(I) 


06(T) 


= Z22(1) *VEL 

= E(1)#03(T) 

PIv2, 

= B(I)*(HwW(1)-B(I) #TAs20) 
= 05(1) *E(1)*.5 


CONTINUE 
RHOG=RHU#GRAVTY 
# #® # # ® SET UP THE FUNCTIONS FUR THE INTEGRALS (PAGE 5 UF NOTES)FUNCT 85 


PIH = 
KPI = 


PI7/2, 
KAR#PI 


EVALUATE INTEGRALS USING TRAP METHOD 


91 


93 


9% 


I = 1 
INDEX 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 


CONTI 
MASS 
QA = 


OMU 

EDMU 
E2DMu 
E30MuU 
BF = 

BMM = 
ZMA = 
ZwWMA 

EMAS 

ZZWMA 
ZWEMA 
Z2WMA 
E2MAZ 
CONTI 
IF ( 

INDEX 
I2zI 
GO TO 


= } 

TRAP (MACINDEX) »DIFF (1) »KTT (1) » TMASS) 
TRAP (QPART (INDEX) sOIFF (I) sKTT (1) 9QAl) 
TRAP (C1 (INDEX) oOIFF (1) oKTT (1) CEA) 
TRAP (C2 (INDEX) »DIFF (1) oKTT (1) »CE2A) 
TRAP (IPART (INDEX) sDIFF (1) sKTT(1I) 9 IAA) 
TRAP (D1 (INDEX) »OIFF (1) 9K TT (1) 9OMUA) 
TRAP (D2 (INDEX) »DOIFF (1) »KTT (1) 9EDMUA) 
TRAP (N3 (INDEX) »DIFF (1) »KTT( 1) 9E2DMUA) 
TRAP (N& (INDEX) oDIFF (1) »9KTT (1) s9E3DMUA) 
TRAP (DS (INDEX) »DIFF (1) oKTT(I) 9 BFA) 
TRAP (N6 (INDEX) oDIFF (1) »KTT (1) »BMMA) 
TRAP (Z1 (INDEX) oDIFF (1) sKTT (1) 9 ZMAA) 
TRAP (72 (INDEX) oDIFF (1) oKTT (1) 9 ZWMAA) 
TRAP (73 (INDEX) »OIFF (1) oKTT(1) 9EMASA) 
TRAP (74( INDEX) »DIFF (1) »KTT(1) 9 ZZWMAA) 
TRAP (75 (INDEX) »DIFF (1) »9KTT (1) 9 ZWEMAA) 
TRAP (76( INDEX) »DIFF (1) »KTT(1) »9Z2WwMAA) 
TRAP (77 (INDEX) »DIFF (I) »KTT (I) s9ECMAZA) 


NUE 

= MASS ¢« TMASS 

QA + QAl 

IA + TAA 

CE + CEA 

cE2 + CE2A 

OMU + DMUA 

= EDM! + EDMUA 

= E20MU «+ E20MUA 
= E30MU «+ E3DMUA 
BF + CHOG*BFA 

BMM + RHOG#®BMMA 
ZMA+ZMAA 

= ZwMA+ZWMAA 

= EMAS+EMASA 
ZZwMA+ZZWMAA 
ZWEMA+Z2WEMAA 
Z2wMA+Z2WMAA 
E2MAZ+E2MAZA 


NUE 

1.GE.11)GO 10 92 
2 INDEX+KTT (I) =} 
ol 

91 


92 CONTINUE 
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FUNCT 77 
FUNCT 78 
FUNCT 79 
FUNCT 80 
FUNCT 81 
FUNCT 82 
FUNCT 83 
FUNCT 84 


FUNCT 86 
FUNCT 87 
FUNCT 88 
FUNCT 89 
FUNCT 90 
FUNCT 91 
FUNCT 92 
FUNCT 93 
FUNCT 94 
FUNCT 95 
FUNCT 96 
FUNCT 97 
FUNCT 98 
FUNCT 99 
FUNCT100 
FUNCT101 
FUNCT102 
FUNCT103 
FUNCT104 
FUNCT105 
FUNCT106 
FUNCT107 
FUNCT108 
FUNCT109 
FUNCT110 
FUNCT111 
FUNCT11l2 
FUNCT113 
FUNCT114 
FUNCT115 
FUNCT116 
FUNCT117 
FUNCT118 
FUNCT119 
FUNCT120 
FUNCT121 
FUNCT122 
FUNCT123 
FUNCT124 
FUNCT12S 
FUNCT126 
FUNCT127 
FUNCT128 
FUNCT129 
FUNCT130 
FUNCT131 
FUNCT132 
FUNCT133 
FUNCT134 
FUNCT135 


Cc # * @ #@ # # # CALL COMPUT TO FIND THE VALUE OF NL AND FL USING FUNCT136 


Cc THE VALUES OF THE ABOVE INTEGRALS FUNCT137 
CALL COMPUT (xX) FUNCT138 

Cc FUNCT139 
IF (NPRINT.LT.3) GO TO 11} FUNCT140 

IF (NPRINT.FEQe3) GO TO 108 FUNCT141 

IF (NPRINT.EQ.4)GO TO 108 FUNCT142 
WRITE (6997) (IPART(I) 9121 oNUM) FUNCT143 
WRITE (6998) (QPART(I) »I21 9NUM) FUNCT144 
WRITE(6999) (C1(I) 9 IT=1 9NUM) FUNCT145 
WRITE (69100) (C2(1I) »I=l NUM) FUNCT146 
WRITE (69161) (C3(I) »IT=2 oNUM) FUNCT147 
WRITE (69102) (Ol (1) 9 I=] 9NUM) FUNCT148 
WRITE (69103) (D2 (1) »J=1 NUM) FUNCT149 
WRITE (69104) (03(1) 9 I=19NUM) FUNCT150 
WRITE(6910S5) (D04(1) o2=) NUM) FUNCT151 
WRITE (69106) (05(1) 5 I=] NUM) FUNCT152 
WRITE(69112) (D6(1) 9 I=1 9NUM) FUNCT153 

WRI FE (69113) (21 (1) » T=1 9NUM) FUNCT154 
WRITE (69114) (Z2(1) 9 T=] 9NUM) FUNCT1S55 
WRITE (60115) (23 (1) 9 L=1 »NUM) FUNCT156 
WRITE (69116) (24(1) 9 T=1 9NUM) FUNCT157 

WRI FE (69118) (25 (1) 9 T= oNUM) FUNCT158 
WRITE (69119) (26 (1) 9 T=] 9NUM) FUNCT1S9 
WRITE (69120) (27 (1) » T= 9NUM) FUNCT160 
WRITE (69107) KPI »RHOG,PIH FUNCT161 

108 WRITE (69109) MASSsCINT»QA9CE,CE2,CE3 FUNCT162 
WRITE (69121)1A FUNCT163 

121 FORMAT(* IA #9E10.4) FUNCT164 
WRITE (69110) OMUsEDMU» E20MU 9 E 3UMU 9 BF 9 BMM FUNCT165 
WRITE (69117) ZMAyZWMA sEMAS 9 ZZWMA, ZWEMA 9 ZZWMA gE2MAZ FUNCT166 

Cc ## # @ &@ & @ © @ © FORMATS # #@ #e #@ toe ee ee FUNCT167 
96 FORMAT (" CPART(1)%510(2XE104)) FUNCT168 
97 FORMAT(" IPART (I) %910(2x9E1004)) FUNCT169 
98 FORMAT(" OPART (I) %510(2X0E100¢4)) FUNCT170 
99 FORMAT(" Cl %910(2X%9E1004)) FUNCT171 
1@0 FORMAT(" C2 %910(2X9E10.4)) FUNCT172 
101 FORMAT(" C3 %510(2XsE10.4)) FUNCT173 
102 FORMAT(" Ol %510(2X2E1004)) FUNCT174 
103 FORMAT(" p02 %o10(2X9E10¢4) ) FUNCT175 
104 FORMAT(" D3 5 10(2X9E1004) ) FUNCT176 
105 FORMAT(" 04 %910(2X%9E1004)) FUNCT177 
106 FORMAT(" DS %510(2X9E1004)) FUNCT178 
122 FORMAT(" 06 %910(2x%9E1004)) FUNCT179 
107 FORMAT(" KPHI "sE100%9"RHUG "E1004" PHIH 961004) FUNCT180 
109 FORMAT(" MASS "“sE10c49" CINT "“sEL0c4%9" QA “sE10c49" CE "sE100%s FUNCT181 
~ @0CE2 5610.40" CE3 9E10.4) FUNCT182 
110 FORMAT (" DMU "9610.45 EDMU % 9610049" E2DMU 9610049" E3DMU ty FUNCT183 
#E10.45" BF 5E10.45" BMM "5E10.4) FUNCT184 
113 FORMAT(SH 21 910(2X5E10.4)) FUNCT185 
114 FORMAT(GH 22 910(2X%sE1004)) FUNCT186 
115 FORMAT(4H Z3 910(2X%sE1004)) FUNCT187 
116 FORMAT(4H 24 910(2X5E10-4)) FUNCT188 
118 FORMAT(4H 25 910(2X5,E10.4)) FUNCT189 
119 FORMAT(4H 26 o10(2XxX5E10.4)) FUNCT190 
120 FORMAT (4H 27 910(2X5E10.4)) FUNCT191 
117 FORMAT(SH 2MA 9€10.45,6H ZWMA 9610.04 96H EMAS 961004 FUNCT192 
es TH Z2WMA 5E10.4 97H ZWEMA »5E10.497H ZZWMA 9E10.4, FUNCT193 

S 7H E2MAZ 5£10.4) FUNCT194 
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lil 


20 


CONTINUE 

RETURN 

END 

SUBROUTINE COMPUT (xX) 

OIMENSION x (6) 

REAL KAR KPI 

REAL NLoMASSoNCGoMoIToITAagK oMAoMMAX oN 
INTEGER END 


COMMON /SHIP/ MASS sCINTsGAsCE sCE2CE3, DMU »EDMU,E2DMU 9 E 3DMU 9 BF » BMM, COMPUT 


NL oFLoIAsE (120) 


~ COMMON JCONST/ NCGeECGoPI sDPReRPD »GRAVTY oRHO 9K »NUM MA (120) oCDo TAs 


ry B(120) sBETAgHW(120) 5 TZ 9>DRAGoWoXDoToXPoMolTo 
¢ DELTAS + TXoEST (120) 6CoROsKARyMMAX (1 0) sTEST(120) » 
N (120) »PHALF 

*COMMON/OUT /NPRINTs NPLOT »END 

COMMON /TEPMS/ TloT20T30T4%oTS9T69T79T8 

COMMON /WAVE/ R(120) »PT (120) » ZMAy ZWMAgEMAS 5 ZZWMA 9 ZWEMA y ZOWMA g 

E2MAZ»ZWDOT (120) 
* COMMON JTEST/ VMA 


Cx6 = COUS(x(6)) 
SxX6 = SIN(X(6)) 
wO = KeC 


CONS] = RO#WO#WOFCX6 

CONS2 = (KPT#RHO#PTH/TA) /CX6 

CONS3 = RO#WOFKECX6eSKE 

CONSS = RO#WOFK#®CX6#CX6 

TERMI = X(1) #®CX6 

TERM2 = X(2)#SX6 

UVNUM = (X(1) #CX6=(X (2) —ZWDOT (NUM) ) #SX6) # 

© (X (1) ®SX6—X (3) #E (NUM) + (X (2) -ZWDOT (NUM) ) #CX6) 


ZMA = ZMA#K (3) #SK6 
ZZWMA = ZZWMA*X (3) #SX6 
ZWMA = ZwWMA#®CUNS1 

EMAS = EMAS#CONS] 

DMU = DMU*CONS2 

EDMU = EDMU#CONS2 

CE = CE*CD#RHO 

CE2 = CE2#cD#RHO 


E2DMU = E2NMU*CONS3 

E30MU = E3DMU*CONS3 

ZWEMA = ZWEMA®CUNS4 

Z2WMA = Z2WMA*CONS4 

Tl = QA#X (3) #(TERMI-TERM2) 
Tl = Tl + 2ZWMA = EMAS 

T2 = EDMU 

T3 = CE2 

T4 = MA(NUM) ®E (NUM) ®UVNUM + E2MAZ + E3DMU - Z2WMA + BMM 
NL = Tl + T2 + T3 + TG + BMM 
TS = MASS®X (3) # (TERM2=TERM1) 
T5 = 15 + 2WMA = ZMA 

T6 = -DMU 

17 = -CE 
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FUNCT19S 
FUNCT196 
FUNCT197 
COMPUT 
COMPUT 
COMPUT 
COMPUT 
COMPUT 
COMPUT 


oorouft wh 


COMPUT 
COMPUT10 
ComMPUTI1 
COMPUT 12 
COMPUT13 
COMPUT 14 
COMPUT15S 
COMPUT16 
COMPUT17 
COMPUT18 
COMPUT19 
COMPUT20 
ComPuT21 
COMPUT22 
COMPUT23 
COMPUT 24 
COMPUT2S 
COMPUT26 
COMPUT27 
COMPUT28 
COMPUT29 
COMPUT30 
COMPUT31 
COMPUT32 
COMPUT33 
COMPUT 34 
COMPUT3S 
COMPUT 36 
COMPUT37 
COMPUT38 
COMPUT39 
COMPUT4O 
COMPUT41 
COMPUT42 
COMPUT43 
COMPUT4&4 
COMPUT4&S 
COMPUT46 
COMPUT47 
COMPUT48 
COMPUT49 
COMPUTSO 
COMPUTS1 
COMPUTS2 
COMPUTS3 
COMPUTS4S 
COMPUTSS 
COMPUTS6 
COMPUTS7 


T8 = =-MA(NUM) #UVNUM = E2DMU + ZWEMA 
BF = BF/CX4 


FL2=TS+¢T6+T7+T8-BF 


IF (NPRINT.LT.3)GO TO 30 
25 CONTINUE 
WRITE (6910) NLoFL 
10 FORMAT(" NL = ",E12.60" FL = 9612.6) 
30 RETURN 
END 
SUBROUTINE INPUT 
C# # © @ © @ © DEFINITION OF INPUT VARIABLES 
XA = INITIAL TIME 
KE = FINAL TIME 
HMIN = MINIMUM STEP SIZE 
HMAX = MAXIMUM STEP SIZE 
EPSE = RELATIVE ERROR CRITERIUN USED FOR VALUES OF Y GT A 
EPs 2 ERROR CRITERION IN KUTMER 
A = ABSOLUTE ERROR CRITERIA USED IN KUTMER 
NPRINT = 1 FINAL PRINTOUT 
= 2 MATRIX INVERSE MATRIXoF COLUMN MATRIX»sAND KUTMER 
RESULTS 
3 INTEGRAL VALUES 
4 CALCULATED VALUES-CONSTANT FOR GIVEN INPUT VALUES 


—S it il 


NPLOT = NO PLOT 
= PRINTER PLOT 
END = NUMBER OF RUNS 


M = MASS OF CRAFT 

WwW 2 WEIGHT OF CRAFT 

Z = THRUST COMPONENT IN Z UIRECTION 

= THRUST COMPONENT IN x DIRECTION 

XECG = DISTANCE FROM CG TO CENTER OF PRESSURE FOR NORMAL FORCE 
xP = MOMENT ARM OF PROPELLER THRUST 
xD = DISTANCE FROM CG TO CENTER UF PRESSURE FOR DRAG FURCE 

KA(I)= ADDED MASS COEFFICIENT 
AN ARRAY GIVEN THE VALUE KAR wHICH IS READ IN 

BM(I) = BEAM AT FREE SURFACE UR AT CHINE 


DRAG = FRICTION DRAG 
K = WAVE NUMBER 
RO = WAVE HEIGHT 
NU = WAVE SLOPE 
NUM = NUMBER OF STATIONS 
BL = BOAT LENGTH 


LAMBDA = WAVE LENGTH 
RG = RADIUS OF GENERATION IN FEET 
T = PROPELLED THRUST IN LBS 
GAMMA = PROPELLER THRUST ANGLE IN VEGREES 
DELTAS=STATION SPACING IN FEET 
ECG = LONGITUDINAL CENTER OF GRAVITY 
NCG = VE?TICAL CG 
BETA(I) = DEAD RISE 
NO(I) = HEIGHT OF MEAN BUTTOCK 
RHO = DENSITY OF WATER 
GRAVTY = GoAVITY FT/SEC##2 
DPR = DEGPEES PER RADIAN 
RPD = RADIANS PER DEGREE 
3.14159 2 oo 0 o © 


DADMDAANAANAAANANANAANAMDDMAMNMNAMNANNANANNANANANANNANNANANAANAANANAANAANANANA 
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COMPUT58 
COMPUTS9 
COMPUT60 
COMPUT61 
COMPUT62 
COMPUT63 
COMPUT64 
COMPUT65 
COMPUT66 
COMPUT67 
COMPUT68 
INPUT 2 
INPUT 3 
INPUT 4 
INPUT 5 
INPUT 6 
INPUT 7 
INPUT 8 
INPUT 9 
INPUT 10 
INPUT 11 
INPUT 12 
INPUT 13 
INPUT 14 
INPUT 15 
INPUT 16 
INPUT 17 
INPUT 18 
INPUT 19 
INPUT 20 
INPUT 21 
INPUT 22 
INPUT 23 
INPUT 24 
INPUT 25 
INPUT 26 
INPUT 27 
INPUT 28 
INPUT 29 
INPUT 30 
INPUT 31 
INPUT 32 
INPUT 33 
INPUT 34 
INPUT 35 
INPUT 36 
INPUT 37 
INPUT 38 
INPUT 39 
INPUT 40 
INPUT 41 
INPUT 42 
INPUT 43 
INPUT 44 
INPUT 45 
INPUT 46 
INPUT 47 
INPUT 48 
INPUT 49 


ANNNANAANANAANN 


a0 


aa0 


EST(I) = STATION POSITION 


START = 
RISE = 


eo * @ #@ © # # » IC OPTIONS 


REAL IT 9K LAMBDA. MoMA o>MMAX gNUsNoNCGgNO MASS oNL og LAgKAR 


COMMON /CONST/ NCGeECGoPI sDPRaRPD sGRAVTY »sRHO 9K »NUM9MA(120) »CDoTAg 
B(120) sBETAsHW (120) > TZyDRAGgWoXDoToXPoMolTo 
DELTAS 9 TX9 EST (120) »CoROoKAR 9 MMAX (170) » TEST (120) 9 


& 
6 


C 
COMMON /IN/ BM(120) 981 (120) »VELIN 
COMMON /IN2/ NO(120) »XA9XE9HMAXsHMINSA (6) o0EPSE (6) sLAMBDA 


C) 


eeamgeenonedee 


IC(1l) =1 USE WAVE Z DISTANCE IN COMPUTING LIFT COMPONENT 


OF NL AND FL 


INTEGER END 


N(120) »PHALF 


* COMMON JSHIP/ MASS oCINT 9QAoCE 9 CE2 0CE3 yDMU»EDMU 9 E2DMU 9 E3DMU 9 BF 9 BMM y 


NL oFL eo ITAgE (120) 


COMMON/OUT /NPRINT »NPLOT 9 END 


START TIME OF THE RAMP FUNCTION FOR SEA WAVE 
DURATION OF THE RISE FROM ZERO TO ONE OF THE RAMP 


COMMON /ACCEL/ XACCL »BWACL »CGACL BL 


NAMELIST/HSP/A gNPRINT »NPLUT gEND pWo Kh 9 TZ9 TX oXECGoXP oXDy 
DRAG »RGo Ts GAMMA sECU »NCG 9 KAR» ROgLAMBDA pNUM,BETASEST 
9X Ao XE »9HMINgHMAX EPS o VELIN 


DATA A bie 015200015.00001501,,.V001».00001/7 


DATA NPRINT»NPLOTsEND/1 9191/7 


DATA WoBLoTZoTX9XECGyXP9XD 9DRAG»RU»sLAMBDA 9RGo T 9 GAMMA 9 
ECGeNCGoKAR /166 930759640009 00416922659 095629280 00% 


2232590.091,0/ 
“DATA NUM »BETASEST 477520.05 


0.00905 .031255 6062509 0093755 2125005 0156255 6187509621875, 

025000 5.281259 631250 9 6343759 6375005 0406255 6437509 0468755 

0500995531255 656250 9 0593759 6625005 6656255 2687596718759 

0750005. 781255 681250 5 0843759687500» 690625 937509 69687591.000 

1.06250 .1.1250051.1875091 62500091 .312551.3750091.4375, 

1.50051 ,5625 91 662591 668759167591) 0812551087591 09375 52005 

200625 5 2.125 926187592 025920312992 037592643715 9265924562592 06255 

2 68759267592 081259208150 9209315930059 30062553012593,1875, 

30250053.3125 9363759304375 030593056255 3.6255 30687543615 / 
DATA XAsXE ,HMINsHMAXsEPS /0,0920.090025y0l0015/ 


DATA VELIN /19,62/ 


# @# 9 # # # READ IN AND WRITE OUT KUTMER PARAMETERS AND PROGRAM 


OPTIONS 

READ (S,HSP) 
WRITE (69HSP) 
OO 10 I=194 
EPSE(I) = EPS 


# # # # # # SET UP CONSTANTS 
PI = 3.141592653589 
GRAVTY=32.18 

DPR=57 29577951308 

RPD=.01 7453292519 
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INPUT 50 
INPUT Sl 
INPUT S2 
INPUT 53 
INPUT 54% 
INPUT 55 
INPUT 56 
INPUT 57 
INPUT 58 
INPUT 59 
INPUT 60 
INPUT 61 
INPUT 62 
INPUT 63 
INPUT 64 
INPUT 65 
INPUT 66 
INPUT 67 
INPUT 68 
INPUT 69 
INPUT 70 
INPUT 71 
INPUT 72 
INPUT 73 
INPUT 74 
INPUT 75 
INPUT 76 
INPUT 77 
INPUT 78 
INPUT 79 
INPUT 80 
INPUT 81 
INPUT 82 
INPUT 83 
INPUT 84 
INPUT 85 
INPUT 86 
INPUT 87 
INPUT 88 
INPUT 89 
INPUT 90 
INPUT 91 
INPUT 92 
INPUT 93 
INPUT 94 
INPUT 95 
INPUT 96 
INPUT 97 
INPUT 98 
INPUT 99 
INPUT100 
INPUT1O1L 
INPUT1O02 
INPUT103 
INPUT104 
INPUT10S 
INPUT106 
INPUT107 
INPUT108 


aa0 


IF (EST (NUM) .LT23075) STOP 3 
COMPUTE NO AND BM ARRAYS 


DO 32 I=1,NUM 
IF(EST(I) .GE.0.75) GO TO 30 

NO (T) =-0046875# (100-SQRT (EST (1) /00375-(EST (1) 0075) ##2,0) ) 
BM (1)=.375#SQRT (1.0-(EST(1) 7. 75-10) ##200) 

GO TO 32 


30 NO(I)=0.0 


BM(I) = 0.375 


32 CONTINUE 


Ceeeteee2COMPUTE CONSTANTS AND INITIALIZE ARRAYS 


ANANNQANAANANANANANAN 


oO 


Cc 


IN 


M=W/GRAVTY 

RHO=1.99 

IT=M#RG*RG 

K = 2.%PT/LAMBDA 
C=SQRT (GRAVTY/K) 
NU=RO#K 

PHALF = (PI/2.)#RHO 


BETA = BETA#RPD 

CO = COS(BFTA) 

TA = TAN(BETA) 

DO 60 I=1.NUM 

E(I) = ECG-EST(I) 

N(I) = NCG*NO(I) 

MMAX(1) = KAR®PHALF®BM(I) #BM(1) 
TEST(I) = (20%BM(1)#TA)/PI1 


CONTINUE 
END=END*1 
RETURN 
END 
SUBROUTINE PLUTER(FX»XA»HMAX »LAMBUA, IB yNWAVE 9 IPT) 
PUT: 
FX A TwO DIMENSIONAL ARRAY CONTAINING PITCH AND 
HEAVE VALUES AT EACH TIME STEP 
KA INITIAL TIME 
HMA x TIME INTERVAL» PTIME*HMAX = INTERVAL BETWEEN 
FX VALUES 


LAMBDA WAVELENGTH USED IN CALCULATING PITCH AND 
HEAVE RATIOES 

IB NUMBER OF FX VALUES 

NWAVE START OF VALUES AFTER WAVE IS COMPLETELY ON 


REAL ITsK »LAMBDA 9MoMA oMMAX gNoNCG 
INTEGER END 


DIMENSION FX (29400) »FMIN(2) »FMAX (2) sNVAR (2) 


COMMON /CONST/ NCGsECGoPI sDPR»RPD»GRAVTY »RHO 9K »NUMo9MA (120) sCDo TAs 


e B(120) sBETAsHW(120) »TZsDRAGoWoXDoToXPoMoITs 


DELTAS 9 TX oEST (120) sCoROoKAgMMAX (120) sTEST (120) 5 


N(120) »>PHALF 


@ 
COMMON/OUT /NPRINT sNPLOT »END 


C# # # # ¢ # # # SET UP VALUES FOR PLOT AND CREATE PLOT 
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INPUT109 
INPUT110 
INPUT1L11 
INPUT112 
INPUT113 
INPUT114 
INPUT115 
INPUT116 
INPUT117 
INPUT118 
INPUT119 
INPUT120 
INPUT121 
INPUT122 
INPUT123 
INPUT124 
INPUT125 
INPUT126 
INPUT127 
INPUT128 
INPUT129 
INPUT130 
INPUT131 
INPUT132 
INPUT133 
INPUT134 
INPUT135 
INPUT136 
INPUT137 
INPUT138 
INPUT139 
INPUT140 
INPUT141 
PLOTER 
PLOTER 
PLOTER 
PLOTER 
PLOTER 
PLOTER 
PLOTER 
PLOTER 
PLOTERIO 
PLOTER1L1 
PLOTER12 
PLOTER13 
PLOTERI14 
PLOTERIS 
PLOTERI6 
PLOTERI7 
PLOTERI8 
PLOTERIO 
PLOTER20 
PLOTER21 
PLOTER22 
PLOTER23 
PLOTER24 
PLOTER2S 
PLOTER26 
PLOTER27 


ODINONL WHR 


Cc 


AaAaqanNaanan 


e # 


a 


200 


7e0 


Ao0 


NF UN=2 

* # # # # # SET UP MIN AND MAX LIMITS FOR PLOT 
FMIN(L)=FX(Lol) 

FMIN(2)=FX(201]) 

FMAX(1)=FX(191) 

FMAX(2)=FK(201) 


# ® # # # » SET UP MIN AND MAX LIMIMTS FOR HITCH AND HEAVE RATIO 


FMNP=F xX (2 oNWAVE) 
FMXP=F X (2 »NWAVE) 
FMNH=F x (1 oNWAVE) 
FMXH=F X (1 »NWAVE) 


00 200 I=1,18 

IF CFX(LoT) LToFMIN(L) DFMIN(L)=FX (Lol) 

IF (FX(L oT) .GToFMAX(1) DFMAX(L) ZF X(1 ST) 

IF (FX(29T) LT oFMIN(2) DFMIN(2) =Fx(2oI) 

IF (FX (201) GT eFMAX (2) )FMAA (2) =F X (291) 

IF (I1-LEeNWAVE)GO TO 200 

IF (FX (oT) .LTeFMNH) FMNH3F X (151) 

IF (FX (191) GT eFMXH) FMXH2F X (151) 

IF (FX (201) .LToF MNP) FMNP2F X (291) 

IF (FX (291) GT eFMXP) FMXP2FX (201) 

CONTINUE 

IF (IPT.EQ.9) GO TO 800 

® # # # # COMPUTE RATIOES 

COL3 = (FMXH=FMNH) 7 (2.#RO) 

COLS% = (FMXP=FMNP) /((4,.#PI#RO) /LAMBDA) 

WRITE (49700) COL3,COL4 

FORMAT(LH1L," HEAVE AMPLITUDE/WAVENEIGHT = "sE12.69/92Xo 
© % PITCH AMPLITUDE (2.#PI®wAVEHEIGHT/LAMBDA) = "sE12.6) 


CONTINUE 

NVAR(1)=10H HEAVE 

NVAR(2)=10H PITCH 

Nls2 i 

XO02KXA 

OELX = HMAx 

IF (NPLOT..E(),-1) CALL PLOT2(FXsFMINoF MAX »NVAR yg NFUN N15 IB,X0,D0ELX) 
RETURN 

END 

SUBROUTINE TRAP (F »0X»NPTS»oANS) 


INPUT? 


F ARRAY OF FUNCTIONAL VALUES OF THE INTEGRAND 
Ox THE X INTERVAL BETWEEN VALUES 
NPTS THE NUMBER OF VALUES GIVEN 


OUTPUT? 


999 


ANS THE VALUE OF THE INTEGRAL 


DIMENSION F(NPTS) 
ANS=0.0 

IF (NPTS.LT.2)G0 TO 999 

DO 1 I=lsNPTS 

ANSZANS*F (1) 

ANS=DX# (ANS-005#(F (1) ¢F (NPTS) )) 
CONTINUE 

RETURN 

END 

FUNCTION RMP(T»START»RISE) 
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PLOTER28 
PLOTER29 
PLOTER30 
PLOTER31 
PLOTER32 
PLOTER33 
PLOTER34 
PLOTER3S 
PLOTER36 
PLOTER37 
PLOTER38 
PLOTER39 
PLOTER4O 
PLOTER41 
PLOTER42 
PLOTERS3 
PLOTER44S 
PLOTER4S 
PLOTER4S6 
PLOTERS7 
PLOTER4S8 
PLOTER49 
PLOTERSO 
PLOTERS] 
PLOTERS2 
PLOTERS3 
PLOTERS4 
PLOTERSS 
PLOTERS6 
PLOTERS7 
PLOTERSS 
PLOTERS9 
PLOTER60 
PLOTER61 
PLOTER62 
PLOTER63 
PLOTER64 
PLOTER6S 
PLOTER66 
PLOTER67 


TRAP 
TRAP 
TRAP 
TRAP 
TRAP 
TRAP 
TRAP 
TRAP 
TRAP 
TRAP 
TRAP 
TRAP 
TRAP 
TRAP 
TRAP 
TRAP 
TRAP 
TRAP 
RMP 


2 
3 
4 
5 
6 
7 
8 
9 
10 
ll 
12 
13 
14 
15 
16 
17 
18 
19 


2 


Cc 
Cc 
Cc 
c 
Cc 
Cc 


80 
99 


# ® #® #@ # THIS FUNCTION IS USED TO GRADUALLY IMPLIMENT THE WAVE 


T CURRENT TIME 

START TIME TO START RAMP FRUM 0.0 TO 1.0 

RISE THE LENGTH OF THE RISE FROM 0.0 TO 1.0 
H=0.0 


IF (T LTeSTART)GO TO 99 
IF (RISE-EQ.0.0)GO TO 80 
TOP=T=START 

Hel Ai!) 

IF (TOP. LT RISE) H=TOP/RISE 
GO JO 99 

H=1. 

IF (T.EQoSTART)H=0.5 
RMP=H 

RETURN 

END 
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RMP 
RMP 
RMP 
RMP 
RMP 
RMP 
RMP 
RMP 
RMP 
RMP 
RMP 
RMP 
RMP 
RMP 
RMP 
RMP 
RMP 
RMP 


LISTING OF COMPUTER PROGRAM FOR CALCOMP PLOTS 


an00 


PROGRAM PLTHSP (INPUT »QUTPUT » TAPES=INPUT » TAPE6=OUTPUT 9 TAPE? 9 TAPES) 


ITAPE = 7 

CALL CALPLT(ITAPE) 
STOP 

END 


SUBROUTINE CALPLT(ITAPE) 
DIMENSION TIME (4003) »PITCH (4003) sHEAVE (4003) 
* » [BUF (1000) »BWACL (4003) »CGACL (4003) 


. LOGICAL ACCEL 


10 


CAL CUMP PLOT OF PITCH AND HEAVE VERSUS TIME 


IREAD = S 
READ(IREAD.10) XAXISsYAXISP ,VAXISHohT 
FORMAT (8F 10.0) 

ACCEL = .FALSE. 

READ(IREAD.20) IA 


9 FORMAT (110) 


IF (TA.EQ.1) ACCEL = .TRUE. 

IF (ACCEL) SCEAD(IREAD,10) YAXISbH,YAXISC 

CALL READT (TIME sHEAVE PITCH» BWACL sCGACL »NPTS) 
CALL PLOTS (IBUF 9100097) 

CALL PLUT(0.591.09=-3) 

CALL ESCALE (TIME »XAXISoNPTSo1) 

CALL ESCALE (HEAVE sYAXISH»NPTS91) 

CALL ESCALE (PITCHs YAXISP»NPTS91) 

IF (ACCEL) CALL ESCALE (BWACL 5 YAXISB»NPTSo1) 

IF (ACCEL) CALL ESCALE (CGACL 5YAXISCsNPTS91) 


Nl = NPTS#} 
N2 = NPTS¢2 
N3 = NPTS¢3 


CALL EAXIS(0.0,0.0,1SHTIME IN SECUNDS 9-15 9XAX1IS50.09 
TIME (NL) o TIME (N2) o TIME (N3) 9HT) 

“CALL EAXIS(0e090ceO0o13HHEAVE IN FEET 9139 VAXISHS90009 

e HEAVE (11) sHEAVE (N2) sHEAVE (N3) 9 HT) 
TEMP = TIME (N2) 
TIME (N2) = TIME (N2)/TIME(N3) 
HEAVE(N2) = HEAVE (N2) /REAVE (N3) 

CALL LINE (TIME sHEAVE »NPTS919090) 

TIME(N2) = TEMP 

XNEW = XAXTS¢3. 

YNEw = 1.0 

CALL PLUT (xXNEW50.05=3) 

CALL EAAIS(00490.0,1SHTIME IN SECONDS 9=155XAXIS 90.05 

e TIME (N1) o TIME (N2) 9 TIME (N3) 9HT) 

CALL EAXIS(00090.0e13HPITCH IN RAVo 913s YAXISP 590009 

6 PITCH(NI]) sPITCH(N2) sPITCHI(N3) sAT) 

TIME(N2) = TIME (N2) /TIME(N3) 

PITCH(N2) = PITCH(N2) ZPITCH(N3) 

CALL LINE (TIME sPITCHsNPTS 91 9090) 

IF (~NOTeACCEL) GO TO 30 
TIME(N2) = TEMP 

CALL PLOT (XNEW50.09=3) 


CALL EAXIS(9.0900091SHTIME IN StCUNDS9=-159XAXIS900 OsTIMECN1]) » 


e TIME (N2) » TIME (N3) sHT) 


MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
CALP 
CALP 
CALP 
CALP 
CALP 
CALP 
CALP 
CALP 
CALP 
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CALL EAXTS(2.050.0916HBUW ACCELERATION) 1659 YAXISB990.0 98WACL (N1) »CALP 


* BWACL (N2) sBWACL (N3) sAT) 
TIME (N2) = TIME (N2)/TIME(N3) 
-BWACL(N2) = BWACL(N2) /BWACL(N3) 
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CALL LINE (TIME sBWACLoNPTS91 9000) 


TIME (N2) = TEMP 
CALL PLUT (XNEW,0.05=3) 


CALL EAXTS(000000091SHTIME IN SECUNDS9=159XAXIS90.009TIME (NI) » 
° TIME (N2) 9 TIME (N3) HT) 
CALL EAXIS(0.050.091S5HCO ACCELERATION15, YAXISCs90.0,CGACL(N1) » 


* CGACL (N2) »CGACL (N3) sHT) 


30 


20 


10 


19 


TIME (N2) = TIME (N2)/TIME(N3) 
CGACL(N2) = CGACL (N2) /CGACL (N3) 
CALL LINE (TIME »CGACL »NPTS»1 90,0) 
CONTINUE 
CALL PLOT (30.090.0,999) 
RETURN 
END 
SUBROUTINE READT (TIME »HEAVE »PITCHs BWACL »CGACL yNPTS) 
DIMENSION x (6) sHEAVE (1) sPITCH(1) 


o 9 TIME (1) sBWACL (1) »CGACL (1) 
I =0 
CONTINUE 
I = Iel 


READ(9) TIME(T) 9 (X(1) 9 1=496) »BWACL (I) »CGACL (I) 

IF (EOF (9))109015 

CONT INUE 

WRITE (6920) TIME(I) 9 (X (J) 9J=496) sBWACL (1) »CGACL(I) 
FORMAT(1H .6(F7e202X)) 

HEAVE(I) = xX(5) 

PITCH(I) = x(6) 

IF (1-GE-40N0) GO TO 10 

GO TO S 

CONTINUE 

NPTS = l-l 

RETURN 

END 

SUBROUTINE EAXIS (XPAGE » YPAGE 9 IBCD »NCHAR gAXLENg ANGLE oF IRSTV 9 
6 DELTAV»sDELTAUsHT) 
DIMENSION 1BCD(1) 


THIS ROUTINE wORKS LIKE THE CALCOMP AXIS WITH THE 
EXCEPTION THAT THE TICK MARKS ARE NOT NECCESSARILY 
EVERY INCH AND THE HEIGHT OF THE CHARACTERS IS INPUTTED 


CALL PLUT(xPAGE » YPAGE 53) 
ISN = ISIGN(1 »NCHAR) 
ISGN = SIGN(1.»DELTAV) 


AMIN = FIRSTV 
x = XPAGE 
Y = YPAGE 
XNUM = FIRSTV-DELTAV 
N = AXLEN/DELTAU 


IF (N@DELTAUCLT-AXLEN) N=Nel 
AMAX = AMIN*(N#DELTAV) 
NDIG = NDIGIT(AMINsAMAX»DELTAUsND) 
CONTINUE 
TEST = (NDIG#HT) + HT 
IF (TEST»GT.DELTAU) HT=HT/2. 
IF (TEST.GT VELTAU) GO TO 10 
AYN = (1¢5#HT) 
BYN = (((NDIG=2) #HT) /204%.9#HT) 
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N = N*l 
TANG = (90.*ANGLE) /57.22958 
ANG = ANGLE/s57,2958 
ST = SIN(TANG) 
cT = COS(TANG) 
S = STN(ANG) 
Cc = COS(ANG) 
00 30 I=leN 


IF (I.EQ.1) GO TO 20 
X = X*DELTAUSC 
Y = Y*DELTAU#S 
CALL PLOT (X9Yo2) 
IF (I-EQeN) GU TO 20 
XT = Xe(,1#CT#ISN) 
YT = Yo(, 1#ST#ISN) 
CALL PLOT (XTsYT92) 
KN = X*¢AYN#CT#ISN=BYN®C 
YN = YeAYN#ST#ISN=BYN#®S 
XNUM = XNUM*¢DELTAV 
CALL NUMBER (XNe YNoHT »XNUM»9 ANGLE 9 ND) 
CALL PLOT (X9Yo3) 
CONTINUE 
XSP = (( (AXLEN/HT) 72.) —( TABS (NCHAR) /2.)) #HT 
YSP = 3.5#HT 
XT = XPAGE * XSP#C ¢ ISN#YSP#CT 
YT = YPAGE * XSP#S + ISN#YSP#ST 
CALL SYMBOL (XT YT oHT»sIBCD, ANGLE 9 LABS (NCHAR) ) 
RETURN 
END 
FUNCTION NDIGIT (AMIN, AMAX 9 ANUM »ND) 


FINDS THE NUMBER OF DIGITS NECCESSARY TO PRINT 
EVEN INCREMENT OF THE FUNCTION ON THE AXIS 


NDIGIT THE NUMBER OF PLACES IN THE ENTIRE NUMBER 
ND THE NUMBER OF DECIMAL PLACES 
ANUM THE VALUE GIVEN TO EACH INCREMENT ON THE AXIS 


IF (ABS (AMIN) oLT.ABS(AMAX)) GO TO 20 
IF (A8S (AMIN) .EQ.ABS (AMAX) eANDeAMAX.NE.0) GO TO 20 
IF (ABS (AMIN) eGT.ABS(AMAX)) GO TO 10 
AMAX = 1, 
AMIN = -1. 
GO TO 20 
AMAX = ABS(AMIN) 
IF (AMAXeLE.1.) GO TO 50 
NOIV = 19 
I SB i 
IF (AMAX/NDIV.LT.1) GO TO 40 
I = I+] 
NDIV = NDIV#10 


= 
wou 


IF (AMAX#NDIV.GT.1.) GO TO 70 
I = [+l 
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NDIV = NDIV#10 
GO TO 60 
NDIGIT = [+2 
ND = I 
DD = FLUAT(ND) 
X = ANUM*(10*#DD) 
IX = X 
IF (X-FLOAT (1X) eLT220001) GO TU 90 
DD = 00+1 
ND = ND¢l 
NDIGIT = NNIGIT#l 
GO TO 80 
CONTINUE 
RETURN 
END 
SUBROUTINE ESCALE (ARRAY »AXLENsNPTS» INC) 


FINDS THE SCALE TO BE USED ON THE AXIS = 
ARRAY MUST HAS THREE UNUSED POSITIONS 
AROAY (NPTS#1) FIRSTV 
ARE AY (NPTS+2) 


VALUES = NUMBERS) 
DELTAU (THE INCREMENT IN INCHES 
BETWEEN TICK MARKS ) 


AROAY (NPTS+3) 


DIMENSION ARRAY (1) 
AMIN = ARRAY (1) 
AMAX = ARRAY (1) 
ISGN = ISIGN(1>INC) 
INC = ITABS(INC) 
DO 10 Y=] eNPTSoINC 
IF (ACRAY (1) eLTe AMIN) AMIN=ARRAY (1) 
IF (ACRAY (1) eGTeAMAX) AMAX=ARRAY (I) 
CONTINUE 
AUNIT = UNIT (AMIN » AMAX o AXLENoNo ANUM) 
CALL AUJUST (AMIN» AMAX pAUNI T > AXLEN No ANUM) 
ARRAY(NPTS+1) = AMIN 
ARRAY (NPTS+2) = ANUM#ISGN 
IF (ISGN.FQe-1) ARRAY (NPTS#1) = AMAX 
ARRAY (NPTS+3) = AUNIT 


IF (ABS (ANUY) 0EQeAUNIT) ARRAY(NPTS*2) = 1.®ISGN 
IF (ABS(ANUM) .EQeAUNIT) ARRAY(NPTS+3) = 1. 
RETURN 

END 


SUBROUTINE ADJUST (AMIN» AMAX gAUNIT » AXLEN 9 Ng ANUM) 


GIVEN AMIN AND AMAX WHICH ARE VISTINCT VALUES» ADJUST 
THEM SO THAT THEY ARE EVEN MULTIPLES OF AUNIT 


K=1 

MIN = AMIN/ANUM 

IF (AMIN eLT.MIN#ANUM) MIN 
AMIN = MIN#ANUM 

MAX = AMAX/ANUM 

IF (AMAXeGT.MAX#ANUM) MAX 
AMAX = MAX#ANUM 

TERM = AMINe (N=K) #ANUM 

IF (TERMeLT,AMAX) GO TO 20 


MIN-1 


MAX+] 
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K = Kel 

GO TO 10 

20 AUNIT = AXLEN/(N-Ke1) 
N = AXLEN/AUNI To] 
RETURN 


“END 


FUNCTION UNIT (AMIN, AMAX 9 AXLEN oN »ANUM) 


FINDS THE INCREMENT BETWEEN VALUES TO BE USED UN THE 
AXIS IN AS FAR AS LABELING THE TICK MARKS 


FINOS THE NUMBER OF DIVISIONS TO BE MADE ON 


THE AXIS 


FINDS THE SIZE IN INCHES OF THESE DIVISIUNS 


IF (AMINoNE.AMAX) GO TO 10 


AM 
AM 


IN = AMIN-1 
AX = AMAX+l 


10 IF (AMAX.LT.1-AND.AMIN.GT.-1)GU TO 110 


30 


40 


110 


MIN 
MAX 
IF 


= AMIN 
= AMAX 
(AMAX.GToMAX) MAX=MAX¢] 


IF (AMINeLT.MIN) MIN=MIN-1 


IF (M 
IF (M 
N 


IN.LT.0) NwID 
IN.GE.9) NwID 
UM = 10 


MAX=MIN 


IF (NWID.-LT2=NUM) GO TO 60 


NU 
GO 
N = 


M = NUM#10 
To 40 
NWID7 (NUM/10) 


MAXI ABS (MIN) 


IF (MINoL Te eANDeMAX.GTe0) GO TO 70 


IF ( 


N@ (NUM/10) eL TeNwWID) N=Nel] 


ANUM = NUM/10. 
AUNIT = AXLEN/N 


GO TO 
0 NN = 
IF (NN®(NUM/10) .LT.» TABS (MIN) ) NN 


N = 

IF O(N 
N = 

ANUM 
AUNI 
GO T 
NUM= 


160 
TABS (MIN) / (NUM/10) 


MAX (NUM/10) 
#(NUM/10) eLToMAX) N = Nel 
NeNN 
= NUM/10. 
T = AXLEN/N 
O 160 
10 


120 IF (AMAX*NUM.GT.1) GO TO 130 


130 
140 


150 


NUM 
GO T 


= NUM#10 

0 120 

1e/NUM 

AMIN®NUM 

AMAX#NUM 

IF (AMIN@®NUM.LT NI) NI=ENI=1 
IF (AMAX#NUM.GT.N2) N2=N2¢1 
IF (N1.NE.N2) GO TO 150 


AMIN = AMIN-UNITT 
AMAX = AMAX=UNITT 
GU T9 140 

N = N2-NnI 


ANUM = UNITT 


= NNel 


IF (AMIN. LTO. AND. AMAXoLT.0) N=N1=N2 
IF (AMIN LT oOcANDeAMAX eGE VU) N=Ne=N1 


AUNIT = AXLEN/N 


73 


JUST 
JUST 
JUST 
JUST 
JUST 
JUST 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
ONIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 
UNIT 


160 IF (NeGTeS) GO TO 170 
5 N = N#2 
ANUM = ANUM/2. 
AUNIT = AUNIT/2. 


GO TO 169 
170 UNIT = AUNIT 
RETURN 
END 
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